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Abstract. The proposed development is an attempt to en- 
hance aerosol retrieval by emphasizing statistical optimiza- 
tion in inversion of advanced satellite observations. This op- 
timization concept improves retrieval accuracy relying on the 
knowledge of measurement error distribution. Efficient ap- 
plication of such optimization requires pronounced data re- 
dundancy (excess of the measurements number over num- 
ber of unknowns) that is not common in satellite observa- 
tions. The POLDER imager on board the PARASOL micro- 
satellite registers spectral polarimetric characteristics of the 
reflected atmospheric radiation at up to 16 viewing directions 
over each observed pixel. The completeness of such observa- 
tions is notably higher than for most currently operating pas- 
sive satellite aerosol sensors. This provides an opportunity 
for profound utilization of statistical optimization principles 
in satellite data inversion. The proposed retrieval scheme is 
designed as statistically optimized multi-variable fitting of 
all available angular observations obtained by the POLDER 
sensor in the window spectral channels where absorption by 
gas is minimal. The total number of such observations by 
PARASOL always exceeds a hundred over each pixel and 
the statistical optimization concept promises to be efficient 
even if the algorithm retrieves several tens of aerosol param- 
eters. Based on this idea, the proposed algorithm uses a large 
number of unknowns and is aimed at retrieval of extended set 
of parameters affecting measured radiation. 
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The algorithm is designed to retrieve complete aerosol 
properties globally. Over land, the algorithm retrieves the pa- 
rameters of underlying surface simultaneously with aerosol. 
In all situations, the approach is anticipated to achieve a ro- 
bust retrieval of complete aerosol properties including infor- 
mation about aerosol particle sizes, shape, absorption and 
composition (refractive index). In order to achieve reliable 
retrieval from PARASOL observations even over very reflec- 
tive desert surfaces, the algorithm was designed as simultane- 
ous inversion of a large group of pixels within one or several 
images. Such multi-pixel retrieval regime takes advantage of 
known limitations on spatial and temporal variability in both 
aerosol and surface properties. Specifically the variations 
of the retrieved parameters horizontally from pixel-to-pixel 
and/or temporary from day-to-day are enforced to be smooth 
by additional a priori constraints. This concept is expected to 
provide satellite retrieval of higher consistency, because the 
retrieval over each single pixel will be benefiting from coin- 
cident aerosol information from neighboring pixels, as well, 
from the information about surface reflectance (over land) 
obtained in preceding and consequent observations over the 
same pixel. 

The paper provides in depth description of the proposed 
inversion concept, illustrates the algorithm performance by a 
series of numerical tests and presents the examples of prelim- 
inary retrieval results obtained from actual PARASOL ob- 
servations. It should be noted that many aspects of the de- 
scribed algorithm design considerably benefited from expe- 
rience accumulated in the preceding effort on developments 
of currently operating AERONET and PARASOL retrievals, 
as well as several core software components were inherited 
from those earlier algorithms. 
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1 Introduction 

The research presented in this paper aims to develop a 
new retrieval algorithm optimized for deriving maximum 
information content using the data redundancy available 
from advanced satellite observations, such as those from 
POLDER/PARASOL observations. The design of the 
POLDER imager allows collecting rather comprehensive 
characterization of angular distribution of both total and po- 
larized components of solar radiation reflected to space. The 
observations in window channels where the effect of absorp- 
tion by atmospheric gases are minimal are usually used for 
aerosol retrievals. The complete set of such observations 
collected operationally by POLDER/PARASOL over each 
pixel includes angular measurements of both total radiances 
and linear polarization at 0.49, 0.675 and 0.87 pm and angu- 
lar measurements of only total radiances at 0.44, 0.565 and 
1.02 pm. The number of viewing directions is similar for 
all spectral channels and varies from 14 to 16 depending on 
observed geographical location. The completeness of such 
observations is significantly higher comparing to any cur- 
rently operating passive satellite aerosol sensors. In addition 
PARASOL provides nearly global coverage every 2 days. 
Therefore, such complete set of PARASOL observations po- 
tentially provides very valuable basis for enhanced charac- 
terization of global aerosol. 

However, rigorous interpretation of redundant satellite ob- 
servation is a very challenging task. Indeed, the optimized 
inversion requires applying complex multi-variable inversion 
algorithms. Such methods are time-consuming and challeng- 
ing for implementation. This is why the rigorous methods of 
inversion optimization are not generally used for processing 
very large data sets provided by satellite aerosol imagers. In- 
stead, most satellite aerosol retrievals use look-up tables of 
simulated satellite signals pre-computed for some limited se- 
lected scenarios of aerosol and underlying surface combina- 
tions. The modeled scenario that provides the best match of 
the observed radiances is accepted as the retrieved solution. 
With some modification, this strategy is adopted in most of 
aerosol satellite retrievals because it allows rapid operational 
processing of satellite images. For example, it is success- 
fully employed in retrievals of single-view AVHHR (Stowe 
et al., 1997; Mishchenko et al., 1999; Higurashi and Naka- 
jima, 1999), TOMS (Torres et al., 1998), MODIS (Kaufman 
et al., 1997; Tanre et al., 1997; Remer et al., 2005), etc. At the 
same time, applying the same methodology for processing 
observations from imagers with multi-viewing capabilities, 
such as MISR (Diner et al., 1998; Martonchik et al., 1998; 
Kahn et al., 2007, 2009), SEVIRI (Govaerts et al., 2010; 
Wagner et al., 2010; Carrer et al., 2010), or POLDER (Deuze 
et al., 2001; Herman et al., 2005; Tanre et al., 2011), reveals 
some deficiencies of the look-up table retrievals. The multi- 
directional observations have notably higher sensitivity to the 
details of aerosol and surface properties, and the retrieval of 
larger number of parameters is expected. Correspondingly, 


the required comprehensive look-up tables of such obser- 
vations may have larger dimensions and thus be less suit- 
able for operational use. As a result, most look-up table 
based algorithms rely only on the selected sub-sets of the ob- 
servations with highest sensitivity to the aerosol parameters 
and retrieve reduced set of characteristics. For example, the 
current POLDER/PARASOL operational retrieval algorithm 
over ocean (Herman et al., 2005) uses measurements of to- 
tal and polarized reflectances only at two spectral channels 
(0.67 and 0.87 pm) . The look-up table algorithm works ef- 
ficiently for these two channels because they are sensitive to 
the scattering of both fine and coarse mode aerosols. At the 
same time, observations at these two channels are insensitive 
to vertical variability of aerosol and not strongly affected by 
water-leaving radiation. The POLDER/PARASOL retrieval 
over land (Deuze et al., 2001) uses only polarized measure- 
ments of reflected light at the same two channels. Such strat- 
egy is used because the contribution of aerosol into polarized 
reflectance generally dominates over the contribution of the 
land reflectance, while contribution of land surface into total 
reflected radiation is usually comparable or stronger than that 
of aerosol. Therefore, as discussed by Deuze et al. (2001), 
utilization of only polarized observations allows one to de- 
rive aerosol properties and to avoid challenging issue of sep- 
aration of surface and aerosol contributions into the total re- 
flectance. Although this algorithm has successfully provided 
valuable aerosol retrievals from POLDER observations, sev- 
eral shortcomings were identified in the POLDER aerosol 
products. First, PARASOL retrieval over land provides in- 
formation only about fine aerosol particles, because the con- 
tribution of large aerosol (predominantly non-spherical dust) 
over land to the polarized reflectances often is small. In ad- 
dition, the correct interpretation of PARASOL observations 
of desert dust even over ocean surface is challenging due 
to difficulties to model appropriately the light scattering by 
non-spherical particles of desert dust (Gerard et al., 2005). 
Second, since the PARASOL algorithm, both over land and 
ocean, relies on the observations of only two spectral chan- 
nels, the retrieved aerosol spectral properties are not always 
fully consistent with the observations at other channels. 

The retrieval algorithm proposed here fits the complete 
set of PARASOL observation in all spectral channels (with 
the exception of the channels dominated by gaseous absorp- 
tion such as 0.763, 0.765 and 0.910 pm) and including both 
measurements of total radiances and linear polarization (if 
available). Based on this strategy, the algorithm is driven 
by larger number of unknown parameters and aimed on re- 
trieval of an extended set of parameters affecting measured 
radiation. For example, the approach allows the retrieval of 
both the optical properties of aerosol and underlying surface 
from PARASOL observations over land. Also, comparing 
to the current operational PARASOL retrieval, the proposed 
algorithm is designed to provide more detailed information 
about aerosol properties including the particle size distri- 
bution, complex refractive index, parameters characterizing 
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aerosol particle shape and vertical distribution. This setup of 
the aerosol retrieval algorithm is based on accumulated ex- 
perience and current understanding of the high potential of 
using spectral multi-angular polarimetric observations from 
space for improving global aerosol monitoring. Diverse as- 
pects of aerosol retrieval improvements by using advanced 
satellite observations have been already demonstrated and 
outlined in numerous previous studies (e.g. see Kokhanovsky 
and de Leeuw, 2009). For example, the studies by Kahn et 
al. (2007, 2009), Kalashnikova et al. (2005), Kalashnikova 
and Kahn (2006) demonstrated the possibility of deriving not 
only aerosol loading but also some information about aerosol 
particle size, morphology and shape from observations from 
the MISR imager that provides multiple view observations of 
total reflectance in 9 directions in 4 spectral channels (0.44, 
0.55, 0.67 and 0.87 pm). These studies have suggested a high 
importance of using multi-angular observation geometry for 
deriving more detailed aerosol information. However, most 
of the known comparisons of the aerosol parameters derived 
from multi-viewing images (such as MISR and POLDER) 
with the aerosol products of single view satellite sensors do 
not indicate clear advantage of multi-viewing observation for 
aerosol monitoring. For example, Kokhanovsky et al. (2007) 
compared the aerosol retrievals obtained from different satel- 
lite platforms over land with ground-based AERONET ob- 
servations and indicated significant differences in currently 
available satellite products and did not reveal any notable 
advantage of one particular satellite sensor. The study sug- 
gested that retrieval indeterminacies are likely part of the 
observed discrepancies, and their reduction will likely be 
aided by new missions incorporating spectral multi-angular 
polarimeters. Indeed, sensitivity analysis of Mishchenko and 
Travis (1997a,b) related the possibility of potential important 
improvements of satellite aerosol retrievals with use of spec- 
tral multi-angular polarization as well as intensity of reflected 
sunlight. Later studies by Chowdhary et al. (2002, 2005) 
demonstrated the possibility of retrieving the detailed aerosol 
properties from spectral multi-angular Research Scanning 
airborne Polarimeter (RSP) over water. This polarimeter is 
an aircraft-based prototype of the APS instrument projected 
to be part of the future NASA Glory mission (Mishchenko 
et al., 2004, 2007). The analysis of RSP observations over 
land by Waquet et al. (2007, 2009a, b) illustrated the possi- 
bility of reliable aerosol retrievals over reflective land sur- 
faces. It should be noted here that Kokhanovsky et al. (2007) 
did not identify any superiority of POLDER results (in- 
cluded into the comparisons) over other satellite imagers. 
This fact has several probable explanations. First, although 
POLDER sensor has multi-viewing polarimetric capabilities, 
the spectral range of POLDER observation is notably nar- 
rower that spectral coverage of some single viewing sen- 
sors, such as a MODIS. Similar remark is valid for com- 
parisons on multi-viewing MISR satellite instrument. Sec- 
ond, the algorithm processing POLDER observations was 
not designed to take full advantage of positive redundancy of 


spectral multi- angle polarimetric observations. Indeed, the 
POLDER retrieval algorithm by Deuze et al. (2001) oriented 
on rapid operational processing uses polarized reflectance at 
two visible spectral channels. Therefore, the positive po- 
larimetric information from other spectral channels, as well 
as any information from total reflectance observation, was 
not used. Waquet et al. (2007) demonstrated that using po- 
larimetric observations over a wider spectral range is essen- 
tial for aerosol retrieval over land from polarimetry obser- 
vations. Waquet et al. (2007, 2009a, b) followed Deuze et 
al. (2001) approach and used only polarized reflectances. At 
the same time, the Waquet et al. (2007, 2009a, b) algorithm 
was driven by a large number of unknowns and was of sig- 
nificantly higher complexity than the POLDER algorithm. 
An algorithm of such level of complexity has never been ap- 
plied to POLDER/PARASOL observations. In addition, one 
could expect that including total reflectance into such an en- 
hanced retrieval scheme could result in additional improve- 
ments of the aerosol retrieval. Indeed, the spectral angular 
measurements of total reflectance are shown to provide valu- 
able aerosol information even over land surfaces (e.g. Mar- 
tonchik et al., 2004; Liu et al., 2004; Kahn et al., 2005; Diner 
et al., 2005; etc.). In addition, rigorous sensitivity studies 
suggest high importance of using observations of both to- 
tal and polarized reflectances for reliable aerosol retrieval 
(Mishchenko and Travis, 1997a,b; Hasekamp and Landgraf, 
2005b, 2007). That is why the retrieval concept described in 
this paper pursues inversion of both total radiances and lin- 
ear polarization measurments and includes implementations 
of several important algorithm refinements. The realization 
of this concept is expected to result in enhancement of com- 
pleteness and accuracy of POLDER aerosol retrieval. 

The presented algorithm developments essentially rely on 
the available positive research heritage from previous re- 
mote sensing aerosol retrieval developments, in particular 
those from the POLDER and AERONET retrieval activities. 
The general inversion scheme will be designed as multi-term 
LSM fitting by Dubovik and King (2000). Such an inver- 
sion strategy allows for the use of a continuous space of so- 
lutions instead of a limited set of predetermined solutions as 
used in look-up table based algorithms. During more than a 
decade the algorithm developed by Dubovik and King (2000) 
was successfully employed for processing observation of 
AERONET of ground-based sun/sky-radiometers (Holben et 
al., 1998). During this period the algorithm has passed no- 
table evolution and several useful modifications were added 
into the inversion procedure. The modifications of that al- 
gorithm were effectively applied for interpretation of coin- 
cident up and down-looking remote sensing observations. 
For example, Sinyuk et al. (2007) conducted the retrieval of 
both atmospheric aerosol and land surface properties from 
a combination of AERONET data with coincident MISR or 
POLDER satellite observations. Gatebe et al. (2010) im- 
plemented the joint retrieval of detailed properties of multi- 
layered aerosol and underlying surface reflectance from a 
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combination of AERONET data and airborne measurements 
by NASA’s Cloud- Absorption Radiometer. The main details 
of the resulting fine-tuned numerical inversion scheme are 
discussed by Dubovik (2004) and below in Sect. 3 of the 
current paper. The modeling of PARASOL observed re- 
flectances is implemented using approaches and computer 
codes developed previously for accurate radiative transfer 
equation solutions and for modeling aerosol single scattering 
and surface reflectance properties. Specifically, the multiple 
solar light interactions with the atmosphere and underlying 
surface are accounted for using the successive order of scat- 
tering radiative transfer code by Lenoble et al. (2007). This 
approach and actual computer code have been used and re- 
fined in POLDER- 1, POLDER-2 and POLDER/PARASOL 
data analysis. The land surface reflectance of total solar radi- 
ance is approximated by the model of Rahman et al. (1993) 
that has already been successfully used for interpretation of 
MISR (e.g. Martonchik et al., 1998), SEVIRI (e.g. Govaerts 
et al., 2010) and RSP (e.g. Litvinov et al., 2010) observa- 
tions. The reflectance of polarized radiation by land surface 
is approximated by using the models developed by Nadal and 
Breon (1999) and Maignan et al. (2009) from the observa- 
tions by POLDER. These models were also validated by RSP 
observations (Litvinov et al., 201 1). The modeling of aerosol 
single scattering properties is adopted from AERONET de- 
velopments. This scattering model seems rather suitable for 
applying to multi-angular polarimetric observations, since 
this model was demonstrated to accurately reproduce the ob- 
servations by ground-based radiometers that have high sensi- 
tivity to the fine features of angular and spectral aerosol scat- 
tering. For example, the software developed by Dubovik et 
al. (2006) for AERONET allows for very fast simulations of 
scattering by non-spherical aerosol particles. As discussed 
by Herman et al. (2005) and Gerard et al. (2005), the ade- 
quate modeling of scattering by non-spherical aerosol parti- 
cles is critical for analysis of PARASOL observations. 

In addition, as a part of the PARASOL aerosol algorithm 
improvement, a new aspect has been introduced into the con- 
cept of satellite data inversion. Specifically, in order to over- 
come some difficulties related to the limited information of 
the PARASOL observations over a single pixel, the retrieval 
is organized as a simultaneous inversion of a large group of 
pixels within one or several images. For example, derivation 
of aerosol properties over bright land is known to be an ex- 
tremely difficult task. The multi-pixel retrieval regime takes 
advantages from known limitations on spatial and temporal 
variability in both aerosol and surfaces properties. Similar 
ideas have already been used in different forms for improving 
satellite retrievals. For example, Martonchik et al. (1998) de- 
rive the surface reflectance properties from a group of near- 
by MISR pixels in a 16 x 16 km region by relying on the 
similarity of aerosol properties over this area. Govaerts et 
al. (2010) have built the SEVIRI aerosol and surface retrieval 
concept assuming rather limited time variability of the land 
surface reflectance properties. Even more explicitly the idea 


was used in studies by Lyapustin et al. (2008) and Lyapustin 
and Wang (2009), who used the limited time variability of 
land surface for screening cloud-contaminated data or the 
limited space variability of aerosol properties for constrain- 
ing aerosol retrieval from MODIS observations. Tempo- 
ral smoothness constraints on surface reflectance variability 
were used by Quaife and Lewis (2010) in the method for in- 
verting linear bi-directional surface reflectance models from 
MODIS observations. Here, the satellite retrieval is designed 
as a statistically optimized simultaneous fitting of the obser- 
vations over a group of pixels implemented under additional 
inter-pixel constraints. Specifically, the variations of the re- 
trieved parameters horizontally from pixel-to-pixel or tem- 
porary from day-to-day over the same pixel are limited by 
the additional a priori constraints, in a similar manner to how 
it is applied in inverse modeling by Dubovik et al. (2008). 
The inclusion of these additional constraints is expected to 
provide retrieval of higher consistency for aerosol retrievals 
from satellites, because the retrieval over each single pixel 
will be benefiting from coincident aerosol information from 
neighboring pixels, in addition to the information about sur- 
face reflectance (over land) obtained in preceding and conse- 
quent observations over the same pixel. 

It should be noted that this paper focuses on detailed im- 
plementation of core ideas for a new PARASOL retrieval al- 
gorithm, however it does not address many aspects of an op- 
erational implementation of the algorithm. For example, is- 
sues such as cloud-screening, retrieval time requirements and 
other important aspects of algorithm implementation for op- 
erational processing are to be addressed in follow-on studies. 

2 General structure of the algorithm 

The general structure of the algorithm is shown in Fig. 1. In 
order to make the algorithm more flexible it is divided into 
several interacting but rather independent modules. Each 
module has rather particular function. The interactions be- 
tween the modules are minimized for a straightforward ex- 
change of a very limited set of parameters. The “Forward 
Model” and “Numerical Inversion” are the two most com- 
plex and elaborated modules in the developed algorithm. The 
organization of the algorithm by modules enhances the flex- 
ibility in algorithm utilization. For example, the “Numeri- 
cal Inversion” module implements quite universal operations 
that have no particular relation to the physical nature of the 
inverted observations. This module can, in principle, be used 
in any other application not related to atmospheric remote 
sensing. The “Forward Model” module does not have such 
universal applicability as the “Numerical Inversion” module. 
Nonetheless, the “Forward Model” module is developed in 
a quite universal way allowing modeling quite a broad va- 
riety of atmospheric remote sensing problems. As a result 
of such organization of the algorithm, it can equally be ap- 
plied (with minimal changes) to inverting observations from 
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General structure of inversion algorithm 



Fig. 1. The general structure of the retrieval algorithm. 


other satellite sensors or from the ground. In addition, such 
an algorithm structure was helpful in adapting physical mod- 
els and computer routine fragments inherited from previous 
AERONET and POLDER developments. 

The following several Sections of the paper provide a full 
description of the “Forward Model” and “Numerical Inver- 
sion” algorithm modules. A number of optional adjustments 
are suggested for setting both aerosol physical model and re- 
trieval scheme. Although the algorithm is tuned for inverting 
PARASOL observations, some aspects of aerosol parameter- 
ization and inversion implementation (in particular a priori 
constraint settings) can be modified and adjusted for optimiz- 
ing the algorithm performance if it is applied to other remote 
sensing observations. For example, two alternative strate- 
gies are suggested for implementing numerical inversion of 
satellite image observations: conventional pixel-by-pixel in- 
version and a new multi-pixel inversion strategy. According 
to this new multi-pixel approach, the retrieval developed as 
simultaneous inversion of a large group of pixels within one 
or several images. Such a retrieval regime takes advantage of 
known limitations of spatial and temporal variability in both 
aerosol and surface properties. 

3 Forward model of POLDER/PARASOL observations 

The aerosol retrieval algorithm is designed to invert 
the POLDER/PARASOL observations acquired in window 
channels shown in Table 1 , that is: the total radiance in 6 win- 
dow channels: 0.44, 0.49, 0.565, 0.675, 0.87 and 1.02 pm, 
and the linear polarization in 3 of these channels: 0.49, 
0.675 and 0.87 pm, reflected by a ground pixel. In each chan- 
nel, observations of the same pixel are performed nearly si- 
multaneously in up to 16 viewing directions (Deschamps et 


Zenith 



Fig. 2. The illustration of the angular convention used for POLDER 
observation modeling. 


al., 1994). It is assumed that the light observed at the top 
of the atmosphere is only linearly polarized. In the polarized 
channels, besides the total reflected radiance, /, the measure- 
ments provide the Stokes parameters Q and U referred to 
axes perpendicular and parallel to the local meridian plane, 
i.e. Q = / p cos(2o') and U = / p sin(2a) where / p is the po- 
larized component of reflected radiance and a is the angle 
between the meridian plane and the polarization direction. 

Let / = (/, Q, U, V) T and E 0 = (E 0 , 0, 0, 0) r stand, 
respectively, for the Stokes’ vectors of the observed elec- 
tromagnetic radiation and of the incident unpolarized so- 
lar radiation; the subscript ‘T” denotes transposition and 
V is assumed to be negligible. The Stokes’ vector 
/ = (/, Q, U, V) T = I (ft o; <po ; (p\\ k) depends on the 
solar zenith angle #o (Mo = cos (#o)), the observation zenith 
angle & \ (/xi =cos (#i)), the solar and observation azimuth 
angles <p$ and <p \ , and wavelength k. Figure 2 illustrates 
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Table 1. List of measured and retrieved characteristics considered in POLDER/PARASOL algorithm. 


POLDER/PARASOL measurements 

MESUREMENT TYPE: 


^(mo;mi;^o;^i;^) = / (©7; 

- I reflected total radiances 

G(mo;mi;<^o;^i;^) = 2(©j; M) 

- Q component of the Stokes vector 

^(mo;mu <pq\<p\\M) = u i®j'’ x i) 

- U component of the Stokes vector 


OBSERVATION SPECIFICATIONS: 
ANGULAR: 


SPECTRAL: 


I (@j; X i), Q (0y ; X() and U (©y; A;) measured in up to 16 viewing directions, that may 
cover the range of scattering angle 0 from ~ 80° to 180°. 


/(0y; A i) measured in 6 window channels A/ = 0.44, 0.49, 0.565, 0.675, 0.87, and 1.02 pm 
Q (0j-; A,-) and U (0^; A/) measured in 3 window channels A; =0.49, 0.675 and 0.87 pm 


Retrieved characteristics 


AEROSOL PARAMETERS 

C v 

dV(ri)/d\nr 


'-'sph 
n(Xi) 
k(X i) 
ho 


- total volume concentration of aerosol (pm 3 /pm 2 ) 

-0 =1,..., N t ) values of volume size distribution in Nj size bins r/, normalized by C v 

- faction of spherical particles 

-0=1, ..., Nx = 6) the real part of the refractive index at every A/ of the POLDER/PARASOL sensor 

- 0 = 1, •••, Nx = 6) the imaginary part of the refractive index at every A/ of the POLDER/PARASOL sensor 

- mean height of aerosol layer. 

Option: algorithm allows the retrieval of multi-component aerosol. In that case, all above parameters are retrieved 

for each aerosol component 


SURFACE REFLECTION PARAMETERS: 

Rahman et al. (1993) MODEL: 

p 0 (A/) - (i = 1, ..., Nx = 6) first RPV BRDF parameter (characterizes intensity of reflectance) 

a: (A/) - (i = 1, ..., Nx = 6) second RPV BRDF parameter (characterizes anisotropy of reflectance) 

0(X i) - (i = 1, ..., Nx = 6) third RPV BRDF parameter (characterizes forward/backscattering contributions) 

/?o(A/) - (i = 1, ..., Nx = 6) fourth RPV BRDF parameter (characterizes hot spot effect) 

NOTE: ho(X) is retrieved only for the observation conditions: ± 7.5° close to backscattering. In other situations, 

it is fixed as related to ho(X) = Po(X). 

Maignan et al. (2009) MODEL: 

B(Xi) - (i = 1 , ..., Nx - 6) free parameter; 

Option: algorithm allows using alternative surface models: Ross-Li model for BRDF and Nadal-Breon model for BPDF. 


the angular convention adapted for storage and processing 
POLDER data. This convention is used in present algorithm. 

The reflected radiance may be written: 

i (m 0 ;h 1 ;wpi;>-) = (l) 

= L [M scat (0; A) + Mreflec (go! Pi ; <P0‘, (Pi ; A)] E 0 + mult, scat., 

where the terms M scat and M re fl e c correspond to the light re- 
flected as a result of single interaction of incident solar light, 
respectively, with the atmosphere and surface. In Eq. (1) it is 
assumed that polarized light is referred to axes perpendicular 
and parallel to the scattering and reflection planes (here, both 
formed by the solar and viewing directions); and the matrix L 
transforms the Stokes’ vector into the plane of observations 
(details are given in Lenoble et al., 2007). 

Under assumption of plane parallel multi-layered atmo- 
sphere, the single scattering term, M scat , at the top of the 
atmosphere can be expressed as: 


M scat (0; A) — 


( 2 ) 


Mo 


Mo + Ml i= ( 


,T '- 1 (1 - e~ mAxi ) P i (©; A) , 

V 7 4 7 T 


where A T; is the optical thickness of the i - th atmospheric 
layer (i = 1, ..., N numbered from the top to the bottom of 
the atmosphere) and r / is the optical depth of the bottom of 
layer i (i.e. r,- = • Ar^); P; (0; A) and co l 0 denote the 

phase matrix and single scattering albedo of the i - th atmo- 
spheric layer, m = l//xo + l//xi . 

The optical properties r ;, P, (0; A) and co l {) of each atmo- 
spheric layer include the contributions of aerosol (character- 
ized in i - th layer by At ; > a , <a>q . and P^ (0; A)), molecular 
scattering (characterized in i - th layer by Ar/ >mo i, co™f = 1 
and PJ 1101 (0; A)) and atmospheric gases (characterized in 
i - th layer by At ; > gas and 60 ^ = 0). The resulting single 
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scattering albedo co l 0 and phase matrix P; (0; X) of the i-th 
atmospheric layer are: 


co l 0 = 


and 


Mqj A Tp a + A T/ ?mo i 
A a “b A Ti mo i + A T i 


(3) 


gas 


p (0; = ^r,, a P»(0; ^) + A I ,, ml PT 1 (9; ^ (4) 


A T i a + A T/ rri ol + A Ti 


gas 


and the extinction optical thickness r of the atmosphere is 
the sum of the corresponding components: 


7 — ^a H~ T mo i + rg as . (5) 

The properties of atmospheric molecular scattering r mo i and 
pmol ^ are we |] known an( j can be calculated with suf- 
ficient accuracy. The absorption of atmospheric gases r gas 
has rather minor contributions in the POLDER/PARASOL 
window channels and can be accounted for using known cli- 
matologies, as well as using available information from an- 
cillary observations. For example, the present development 
uses the same procedure as used in the operational algorithm 
by Deuze et al. (2001). That procedure corrects the water va- 
por absorption using PARASOL measurements in 0.910 pm 
spectral band. The minor absorption from ozone, NO 2 and 
O 2 are accounted using the climatology data. Thus, the most 
challenging part in modeling single scattering properties of 
the atmosphere is the modeling of aerosol contribution, i.e. 
aerosol extinction r a , single scattering albedo coq and phase 
matrix P a (0; X). These properties depend on aerosol mi- 
crophysics: particle size, shape and composition (refractive 
index). All these characteristics are driven by the parameters 
included in the vector of unknowns and correspondingly they 
are retrieved from the observations. 

The single reflection M re fl ec at the top of atmosphere can 
be calculated as: 

Mreflec (P(f PL <Po; <Pi; k) = ^ e~ mr * R (p 0 ; Pj ; (p 0 ; <pi;X) (6) 

where reflection matrix R (/xo; fi\‘, cpo', <p\\ X) describes the 
surface reflection properties in the plane formed by the so- 
lar and viewing directions. For the ocean surface the re- 
flection R (/xo; /xi ; (po; (p\ ; X) is mainly governed by the 
wind speed at sea level as suggested by the Cox-Munk model 
(Cox and Munk, 1954) employed in the currently operational 
POLDER algorithm (Deuze et al., 2001 ; Herman el al., 2005; 
Tanre et al., 2011). In contrast, the reflection matrix of the 
land surface may differ very strongly from location to loca- 
tion. Therefore, in the present algorithm, the key properties 
of the land surface reflectance are included in the set of un- 
knowns and retrieved from the observations. 

As follows from Eq. (1), once the single scattering terms 
M scat and M re fl ec are defined one needs to account for multi- 
ple interactions of scattered light with atmosphere and sur- 
face. In the present algorithm these interactions are ac- 
counted for by rigorously solving vector radiative transfer 


Forward Model 



Fig. 3. The organization of the forward calculations of atmospheric 
radiance measured from a satellite. 


equation. Thus, the forward model of reflected radiances 
measured by POLDER/PARASOL contains three main com- 
ponents: (i) aerosol single scattering, (ii) surface reflection 
and (iii) solving vector radiative transfer equation for ac- 
counting for multiple scattering. The following parts of this 
section will describe each of these components in detail. 

It should be noted that the forward model for repro- 
ducing POLDER/PARASOL observations is designed by 
means of adapting the atmospheric modeling strategies and 
computer routines developed within previous POLDER and 
AERONET activities. At the same time, several important 
modifications required for optimizing the forward modeling 
performance have been implemented in the present PARA- 
SOL algorithm. Specifically, the models of land surface 
reflectance have been introduced into the radiative transfer 
calculations, the number of aerosol parameters driving the 
model has been reduced, the different regimes of the radia- 
tive transfer calculations have been designed for allowing 
faster but less accurate calculations. These and other forward 
model modification allowed the performance of the devel- 
oped “on-line” inversion procedure to attain the standards re- 
quired for operational processing (achieving sufficient speed 
of computations, etc.). 

Figure 3 shows the data flow within the “Forward model” 
block of the algorithm. Three main complementary efforts 
are involved in modeling the atmospheric radiation held ob- 
served by the POLDER sensor: 

- Modeling of single scattering properties of the atmo- 
spheric aerosol. 
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- Modeling of the surface reflectance properties. 

- Accounting for multiple scattering effects using full ra- 
diative transfer model. 

These aspects are described in detail below. 

3.1 Aerosol single scattering properties 


r max £ max 



(A; k; n; r; s ) 


dN ( e ) dN (r) 
d In e d In r 


d In e d In r 


( 8 ) 


= ^2 K® (A; k\ n; n; s k ) 

i = 1 N e k= 1 , N s 


dN (8k) dV in) 
d In s d In r 


The modeling of the aerosol scattering matrices has been im- 
plemented following the ideas employed in AERONET re- 
trieval algorithm by Dubovik and King (2000) and Dubovik 
et al. (2002b, 2006). 

In order to account for aerosol non-sphericity, the atmo- 
spheric aerosol is modeled as an ensemble of randomly ori- 
ented spheroids. Specifically, AERONET operational re- 
trievals use the concept developed by Dubovik et al. (2006) 
and model the particles for each size bin as a mixture of 
spherical and non-spherical aerosol components. The non- 
spherical component was modeled by ensembles of ran- 
domly oriented spheroids (ellipsoids of revolution). Accord- 
ing to this concept, aerosol particles of non-spherical com- 
ponent have size-independent distribution of shapes and the 
modeling of the total aerosol optical thickness r of non- 
spherical aerosol can be written as the following: 

r (A.) = (7) 


r max £max 

J c\ (A; 1c; n; r; e) 

r xmn £min 



dN (e) dN (r) 
d In e d In r 


d In s d In 


r. 


where c\ (A; k; n; r; s ) denotes the extinction cross-sections 
of spherical particle and randomly oriented spheroid, A - 
wavelength, n and k - real and imaginary parts of the re- 
fractive index, e spheroid axis ratio ( s = a/b , a - axis of 
spheroid rotational symmetry, b - axis perpendicular to the 
axis of spheroid rotational symmetry), r - radius of volume- 
equivalent sphere. Correspondingly the mixture of spheroid 
includes the flattened oblate spheroids (s < 1), spheres (e = 1) 
and elongated prolate spheroids (s > 1). The characteristics 
r and e are used here for describing size and shape of the 
ensemble of spheroids. Analogously to the combination of 
a and b, the combination of r and s allows unique defini- 
tion of the spheroid shape. As discussed by Mishchenko 
et al. (1997), the usage of r and e is convenient for sep- 
arating the effect of particle shape and size in analysis of 
aerosol mixture light scattering. Then the functions 
and d 2\ns denote the number particle size and the number 
particle shape (axis ratio) distributions accordingly. 

For performing fast and accurate calculations of aerosol 
extinction by poly disperse non-spherical aerosols, the size 
and shape integration can be approximated by the double 
sum, e.g. 


where 


Kg Xt (A; k; n; r t ; e k ) = 

lnr,+Alnr lne^+Alng 

r C c®( A; k; n; r; 

J J m 

lnr, — Alnr lne^— Alne 


e) 


(9) 


Ak (e) Bi (r) d In s d In r. 


where v(r) is the volume of particle, d J^ ~ 
the volume particle size distribution, A * (e) and B( (r) are 
the functions providing correspondingly the interpolation of 
shape distribution between the selected points Sk and the in- 
terpolation of size distribution over selected points r,-. In 
studies by Dubovik et al. (2002b, 2006), the coefficients 
Ak ( s ) for integrating over axis ratio were assumed as rectan- 
gular functions Ak (s) = const. For approximating size distri- 
bution between the used size bins r; , the trapezoidal approxi- 
mation was chosen by Dubovik and King (2000). Such kinds 
of interpolation are traditionally applied in aerosol applica- 
tions (e.g. see Twomey, 1977), where the functions Bi (r) 
are defined as isosceles triangles. 

It should be noted that in Eq. (8) the integral over size is 
approximated by a sum using values of volume size distri- 
bution dV(r)/dlnr (in place of the number size distribution 
dN(r)/d\nr ) defined in logarithmically equidistant points r;. 
Utilization of both the volume size distribution and logarithm 
of radius was chosen for the convenience of the algorithm 
implementation. In principle, the particle number distribu- 
tions dN(r)/dlnr or dN(r)/dr could be equally used in Eq. (8) 
(e.g. see King et al., 1978; King, 1982). At the same time, 
the usage of both the volume of the particle (instead of num- 
ber) and logarithmic scale in binning of the size distribution 
helps to optimize the approximation given by Eq. (8). First, 
these choices help to improve the accuracy of this approx- 
imation (a smaller number of points N v provides appropri- 
ate accuracy). Second, under this representation, the ker- 
nels Kg Xt (A; k; n; r Sk) for different points r t are closer 
in values. This is one of the favorable conditions for im- 
plementing inversion. Therefore, volume size distribution 
dV(r)/dlnr is often used as retrieved aerosol characteristic 
in the algorithms applied to invert the optical data of high 
sensitivity to aerosol particle size. For example, a similar 
size distribution representation was used in earlier studies by 
Nakajima et al. (1983, 1996) for retrieving aerosol proper- 
ties from ground-based sky-radiometers. In the AERONET 
retrieval, Dubovik and King (2000) represented volume size 
distribution dV(r)/d\nr by N v = 22 points r;. These points 
are equidistant in logarithmic space and cover the size range 
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from 0.05 to 15 pm. This size range was chosen following the 
sensitivity analysis by Dubovik et al. (2000), which showed 
that the aerosol particles of smaller and larger sizes produce 
negligible contribution to AERONET radiometer observa- 
tions. This range of aerosol particle sizes is slightly wider 
than the one used in earlier studies by Nakajima et al. (1996). 
As discussed later in this paper, the size range of aerosol 
is modified for retrieving aerosol from PARASOL observa- 
tions. 

Using the approximation given by Eqs. (8)-(9), Dubovik 
et al. (2006) developed a numerical tool for fast calculations 
of scattering properties of spheroid mixture. The quadrature 
coefficients K^ xt (X; k; n; rp St) for the extinction, as well 
as for absorption cross-sections and scattering matrices have 
been calculated and stored into the look-up tables for a wide 
range of n and k (13 <n <1.1; 0.0005 < k < 0.5). The cal- 
culations were done for spheroids with axis ratio e ranging 
from ~0.3 (flattened oblate spheroids) to ~3.0 (elongated 
prolate spheroids) and for 41 narrow size bins covering the 
size-parameter range from ~0.012 to ~625. The look-up 
tables were arranged into a software package allowing fast, 
accurate, and flexible modeling of scattering by randomly 
oriented spheroids with different size and shape distribu- 
tions. In addition, Dubovik et al. (2006) used the developed 
software and showed that spheroids can closely reproduce 
single-scattering matrices of mineral dust measured in the 
laboratory by Volten et al. (2001). It was shown that scatter- 
ing matrices have rather limited sensitivity to the minor de- 
tails of axis ratio distribution . Therefore, Dubovik et 
al. (2006) have suggested and demonstrated that AERONET 
retrieval may rely on rather simple assumption that shape 
(axis ratio) distribution in the non- spherical fraction 

of any tropospheric aerosol is the same. Based on this con- 
clusion the aerosol scattering model was set in AERONET 
retrieval as a mixture of spherical and non-spherical frac- 
tions, and d j\ne ) obtained by Dubovik et al. (2006) from 
fitting Volten et al. (2001) observations was employed as 
shape distribution for non-spherical fraction. Based on this 
assumption, the integration over s in Eq. (7) can be done 
once and for all for each size bin, and, therefore, modeling 
of aerosol optical properties (r a (X), coq and P a (0; X)) in 
AERONET retrieval is implemented in a particularly conve- 
nient form. For example, for modeling r a (k) one can write: 

(k) — T S ph (k) + ^nons « (10a) 

= J2 ( C 'pl> K ™ 0-; k; n- n) + (1 - C sph )K”“ (X; k ; n; r,)) 

i= 1, N r 

where 

KjJ (A; k\ n\ n) = (10b) 


In rj+Alnr 

- / 

In rj —A In r 


CgJ (A; k ; n; r) 
v(r) 


Bk (r) d In r, 


i^nons 

^ext 


(X; k\ n\ r,-) = 


(10c) 


1 n r,- + A In r 

= / S ‘ (r) 
In r,-— Alnr 


/ 


cL U; k; n; r; s ) dN(s ) 


v(r) 


d In s 


d In e I d In r, 


where C sp h is the fraction of the spherical particles. Note, 
that while the preparation of the core look-up tables 
K™ t ns (X; k; n; r/) required several years of computations, 
the resulting software allows very fast simulations of scat- 
tering by non-spherical aerosol particles. The calculation 
takes a fraction of a second for any realistic combination of 
aerosol size distribution and complex refractive index. At 
present, it is probably the only approach that can calculate 
scattering matrices for non-spherical particles as part of the 
retrieval without relying on look-up tables of scattering ma- 
trices. Once the forward model in AERONET retrieval was 
updated with non-spherical aerosol light scattering modeling 
capabilities by Dubovik et al. (2006) (as shown in Eq. 10) the 
fraction C sp h was included in the set of retrieved parameters 
along with the concentrations for 22 bins of size distribution. 

It is noteworthy that the spheroid model developed by 
Dubovik et al. (2002b, 2006) appeared to be rather useful 
for AERONET and other aerosol remote sensing applica- 
tions. First, the utilization of this model has significantly im- 
proved the AERONET operational retrieval of aerosol with 
pronounced coarse mode fraction (e.g. see Reid et al., 2003; 
Eck et al., 2005; Dubovik et al., 2006). The same model has 
been shown to reproduce adequately the ground-based po- 
larimetric observations of non-spherical desert dust. Specifi- 
cally, the efficient application of the model to the polarimet- 
ric observations has been done by Dubovik et al. (2006) for 
a case study and Li et al. (2009) for an extended series of ob- 
servations. In addition, it was shown that the spheroid model 
allows qualitative reproduction of the main characteristic fea- 
tures of lidar observations of non-spherical desert dust. For 
example, the increase of extinction-to-backscattering lidar 
ratio and a high depolarization of signal regularly observed 
in lidar observations of desert dust, and traditionally associ- 
ated with aerosol particle non- sphericity, can be adequately 
reproduced using a spheroid-based model (see discussion by 
Dubovik et al., 2006). Cattrall at al. (2005) showed that 
lidar ratios calculated from aerosol properties derived from 
AERONET observations using a spheroid model agree well 
with known lidar observations of desert dust. Furthermore, 
Veselovskii et al. (2010) have used the approach suggested by 
Dubovik et al. (2006) and incorporated the spheroid model 
into the algorithm retrieving aerosol properties from lidar ob- 
servations. That is, probably, one of the first attempts to in- 
terpret quantitatively the sensitivity of the lidar observations 
to particle non- sphericity. The non-spherical coarse aerosol 
models derived from climatologies of AERONET retrievals 
had been successfully incorporated into MODIS and SEVIRI 
satellite retrievals (Levy et al., 2007a,b; Govaerts et al., 2010; 
Wagner et al., 2010). The AERONET retrievals are being 
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Trapezium Approximation 


Approximation by Log-Normals 



Particle Radius, \xm 



Fig. 4. The illustration of modeling aerosol size distribution by five (A r = 5) triangular (left) or log-normal (right) size bins. 


used for making accurate calculations of atmospheric broad- 
band fluxes and aerosol radiative forcing. These calculations 
were shown to agree very reasonably with available coinci- 
dent ground-based flux observations in desert regions (Der- 
imian et al., 2008) and globally (Garcia et al., 2008). De- 
rimian et al. (2008) demonstrated that the neglect of desert 
dust non- sphericity in climatic assessment leads to ~10% 
systematic overestimation of cooling of the atmosphere by 
desert dust aerosol on the top of the atmosphere. 

The retrieval algorithm developed here for POLDER/ 
PARASOL uses the same modeling strategy as described 
above. Correspondingly, a solution is sought in continuous 
space of aerosol size distribution parameters, aerosol parti- 
cle shape and complex refractive indices (see Table 1). How- 
ever, due to differences in information content of AERONET 
and POLDER measurements, the retrieved size distribution 
is represented by a smaller number of bins N v . Instead of 
N r = 22, retrieved by AERONET, here N r is reduced to 16 
and even significantly smaller numbers. In order to assure 
that the aerosol model remains adequate and its accuracy 
is acceptable even if number of aerosol bins is small, the 
performance of Dubovik et al. (2006) software was ana- 
lyzed for situations corresponding to POLDER/PARASOL 
measurements with reduced number of aerosol bins. It 
was found that the accuracy of the calculations remains 
practically unchanged if the aerosol size bins correspond- 
ing to very small and very large particles have been elim- 
inated, i.e. A r = 16 covers size range from r m i n = 0.07 pm 
to r max = 10prn. The contribution of smaller and larger 
particles into POLDER/PARASOL observations is negligi- 
ble. Nonetheless, the need for optimizing the software 
of Dubovik et al. (2006) by using even smaller number 
of aerosol bins N r < 10, was identified. It was found 
that sufficiently accurate modeling of POLDER/PARASOL 


observations can be achieved if the shape of each single bin is 
optimized. Specifically, utilizations of log-normally shaped 
bins provided notable improvements if N r < 10. The con- 
ducted small series of calculations suggested some advan- 
tages of modeling aerosol size distribution as superposition 
of the log-normal functions with fixed parameters: 


dV (r) 
d In r 


E 


cyj 

\/ 2 71 Gi 


exp 




(11a) 


= E « 


dVj (r) 
d In r 


i.e. 

Ci = and (lib) 


dvi (r) 1 (in r - In r v ,/) 2 

d In r V 27r 07 2 af 

For example, Fig. 4 illustrates that a small number of log- 
normal bins retains the realistically smooth shape of atmo- 
spheric aerosol size distribution, while applying triangular 
or trapezium approximation leads to appearance of inade- 
quate features (apparent triangle tops) in the size distribution. 
Thus, the new option allowing the usage of log-normally 
shaped bins was included in Dubovik et al. (2006). 

3.2 Modeling surface reflectance 

The reflective properties of ocean surface are modeled anal- 
ogously to the currently operational POLDER algorithm 
(Deuze et al., 2001; Herman et al., 2005; Tanre et al., 2011). 
The Fresnel’s reflection on the agitated sea surface is taken 
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into account using the Cox and Munk model (Cox and Munk, 
1954). The water leaving radiance is nearly isotropic (Voss 
et al., 2007) and modeling shows that its polarization is neg- 
ligible (e.g. Chami et al., 2001; Chowdhary et al., 2006; Ota 
et al., 2010). This term and the white cap reflection are taken 
into account by Lambertian unpolarized reflectances. The 
whitecap reflectance is driven by the wind speed at sea sur- 
face according to the Koepke (1984) model. The seawater re- 
flectance at short wavelengths is not negligible and depends 
on the properties of oceanic waters. Thus, in present model, 
the wind speed and the magnitude of seawater reflectance at 
short wavelengths need to be known a priori or retrieved to- 
gether with aerosol. 

The modeling of the reflectance by the land sur- 
faces has been adjusted to the needs of newly developed 
POLDER/PARASOL retrieval. The aerosol retrieval algo- 
rithm by Deuze et al. (2001) over land relies only on the 
PARASOL measurements of polarized reflectance and, cor- 
respondingly, it does not consider the detailed directional 
scattering properties of total reflectance by land surface. 
Therefore, the “forward model” was modified to account ad- 
equately for both total and polarized properties of surface re- 
flectance. 

In remote sensing applications the effects of directional- 
ity of land surface reflectance are often accounted for by 
semi-empirical models driven by a small number of internal 
parameters. For example, the Ross-Li model (Ross, 1981; 
Li and Strahler, 1992; Wanner et al., 1995) is employed for 
characterization of directional properties of land surface re- 
flectance derived from MODIS observation (Justice et al., 
1998). The Rahman-Pinty-Verstraete (RPV) model (Rahman 
et al., 1993) is successfully used for the analysis of MISR ob- 
servations by Martonchik et al. (1998) and SEVIRI by Gov- 
aerts et al. (2010) and Wagner et al. (2010). The comparisons 
of the models with satellite (Maignan et al., 2004, 2009) and 
aircraft (Litvinov et al., 2010, 201 1) data showed that, gener- 
ally, both Ross-Li and RPV models are comparably capable 
of reproducing the multi-angle observations of land surfaces. 
Since the RPV model was applied more extensively to inter- 
pretation of multi-directional images (e.g. MISR, SEVIRI), 
it has been retained in the present POLDER/PARASOL al- 
gorithm as a primary formulation for modeling Bidirec- 
tional Reflectance Distribution Function (BRDF). Rahman et 
al. (1993) describe BRDF as: 


prpv Oh; <pi; fe <pi) = 


(12a) 


= Po 


F (g) = 


COS k 1 Pi COS k l &2 


(cos Pi + COS P 2 L 


fr F(g) 
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(12b) 


[l + 0 2 — 2 0 cos (n — at)} 
cos oii = cos Pi cos P 2 — sin Pi sin P 2 cos {ap\ — (P 2 ), (12c) 


G — |\an 2 Pi +tan 2 P 2 — 2 tan Pi tan P 2 cos {<p\ — <^2)J • (12d) 

This model provides bi-directional reflectance as a function 
of four empirical parameters: p Q (X) characterizes intensity 
of reflectance, k (X) characterizes anisotropy of reflectance, 
9 (X) characterizes forward/backscattering contributions in 
total reflectance, ho (X) is a “hot spot” parameter. All these 
parameters are considered as unknowns and included in the 
set of retrieved parameters (see Table 1). 

The available airborne and satellite polarimetric observa- 
tions in the visible and infrared showed that the Bidirectional 
Polarization Distribution Function (BPDF) of land surface 
tends to have rather small values (compared to BPDF) with 
no spectral dependence (e.g., Rondeaux and Herman, 1991; 
Nadal and Breon, 1999; Maignan et al., 2004, 2009; Wa- 
quet et al., 2009a, b; Litvinov et al., 2010, 2011). Most theo- 
retical models developed for approximating observed BPDF 
are based on the Fresnel equations of light reflection from 
the surface. For example, Nadal and Breon (1999) have 
proposed simple two-parameter non-linear function of the 
Fresnel reflection for characterization of atmospheric aerosol 
over land surfaces based on POLDER observations of land 
surface reflectance. Recently, Maignan et al. (2009) have in- 
troduced a new linear BPDF model with only one free pa- 
rameter and demonstrated that this simple model allows a 
similar fit to the POLDER measurements as more complex 
non-linear model by Nadal and Breon (1999). This model 
has been used in the present POLDER/PARASOL retrieval 
algorithm as a primary model of polarized reflectance of land 
surface. 

The model by Maignan et al. (2009) describes the polar- 
ized reflectance as: 

Ps Oh; <PW <P2 ) = PMaignan F U («;, n) , (13a) 

where BPDF is given as a linear function of polarized com- 
ponent F \2 (oii, n ) of the Fresnel reflection matrix (depen- 
dent on incident angle cq- and refractive index n ) multiplied 
by an empirical coefficient: 

. ... , _ B exp (-tan (of,-)) exp (-v) 

PMaignan (h.» Vli tyl) — A ( \ ' O^b) 

4 (Po + Pi) 

The attenuation term exp(— z) reflects the observed tendency 
of decreasing polarized reflectance with increasing vegeta- 
tion cover, where z is the Normalized Difference Vegetation 
Index (NDVI). The NDVI value z was obtained from the 
reflectance measurements concomitant with the polarization 
observations, B is a free parameter that should be chosen to 
fit the observed BPDF angular dependence. 

Thus, in the present algorithm the land surface reflectance 
properties are modeled using Eqs. (12) and (13) for simulat- 
ing BRDF and BPDF accordingly. However, since these for- 
mulations are semi-empirical and derived completely inde- 
pendently, one needs to exclude physically unrealistic com- 
binations of BRDF and BPDF. Therefore, in order to assure 
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that the surface reflectance of polarized radiation in any ge- 
ometry does not exceed the reflectance of total radiation, the 
reflectance matrix R (/xo; /x i; <po; <p\; X) of the land surface 
is represented as a sum of two surface reflection phenomena: 

R = R diff + R spec? (14a) 


or either BRDF or BPDF models. For example, in the 
present algorithm version, the BRDF can be simulated using 
the Ross-Li model and BPDF can be modeled using Nadal- 
Breon (1999) formulation. Nonetheless, the utilizations of 
these secondary BRDF or BPDF models will not be dis- 
cussed in detail. 


where the matrix for diffuse unpolarized reflectance Rdiff is 
modeled using Eq. (12): 


3.3 Full forward radiative model 


R diff (p 0 ; Pi ; (po; <p \ ; k) = prpv 


(\ 000\ 
0000 
0000 
\0000 / 


(14b) 


and the matrix for specular reflectance R spe c is mod- 
eled as matrix of Fresnel reflectance F (a;, n ) scaled by 
PMaignan (#i ; p\\ (pi) empirical coefficient: 


R spec Oh; <PW ^25 VX X) = PMaignanF (Off , n) , (14c) 


where the Fresnel matrix F (cq-, n ) defined for Stokes pa- 
rameters referred to directions parallel and perpendicular to 
the reflection plane can be written as (e.g. see Lenoble et al., 
2007): 


F (off, n) = 


/ r, 2 + r 2 r 2 - r 2 0 0 \ 

r 2 - r 2 r 2 + r 2 0 0 

0 0 2 ri r r 0 

^ 0 0 0 2 ri r r y 


(15a) 


The Fresnel matrix F ( a,- , n ) depends on incident angle d{ 
and refractive index n. The coefficients r r and r\ are defined 
as 


sin (a T — di) tan (a r — cq) 

r r = and r\ = , 

sin ( a v + oii) tan ( d r + di ) 


(15b) 


where the refraction angle d x is related to d\ through the 
Snell-Descartes refraction law given as: 


sin (d() = n sin (a r ). 


(15c) 


The straightforward analysis of the above equations shows 
that the definition of R (/x o; /x i ; (po; <p\; X) given by Eq. (14) 
secures the physically correct ratio between polarized and to- 
tal radiation components (Rij < R\ \ ) for any combination of 
Prpv and PMaignan- 

Thus, in the present algorithm, the BRDF and BPDF prop- 
erties are driven by four free spectral parameters of RPV 
model (p 0 (X), k (X), 9 (X), ho (X)) and one free (generally 
spectrally dependent) parameter B (X). All these parame- 
ters have been added into the retrieved vector of unknowns 
as shown in Table 1 . 

It should be noted, however, that the Rahman et al. (1993) 
and Maignan et al. (2009) formulations, chosen as primary 
models for BRDF and BPDF in the present algorithm, have 
limited accuracy (e.g. see Litvinov et al.,2010, 2011). There- 
fore both the forward model and inversion module of the 
present algorithm have an assumed option of changing both 


Accounting for multiple scattering effects in the atmosphere 
is implemented by the successive order of scattering radiative 
transfer code (Lenoble et al., 2007) that was used in PARA- 
SOL operational retrievals (Deuze et al., 2001 ; Herman el al., 
2005; Tanre et al., 2011). The code provides full information 
about the atmospheric radiation field under the assumption of 
the plane parallel atmosphere. In order to reduce calculation 
time for inverting PARASOL observations, the V component 
in the Stokes vector has been neglected. Hansen (1971) has 
demonstrated that the radiation properties measured by pas- 
sive remote sensing exhibit negligible circular polarization of 
the electromagnetic field. The developed version of succes- 
sive order of scattering radiative transfer code allows calcula- 
tions of atmospheric radiances for several Nk aerosol compo- 
nents. Each aerosol component can be described by defined 
vertical profile of spectral extinction k ex t,£ k) and altitude 
independent phase matrix (0, X) and single scattering 
albedo coq (X). In the present set up of the aerosol retrieval 
code these optical properties are detremined based on mi- 
crophysical model of atmospheric aerosol. Correspondingly, 
only parameters describing aerosol microphysics are directly 
included in the set of retrieved parameters listed in Table 1 . 
Specifically, the vertically invariant P& (0 , X) and coq (X) are 
driven by: the shape of the size distribution dV kin) / hflnr giv- 
ing the aerosol particle volume in the total atmospheric col- 
umn per unit of surface area (in the unites of pm 3 /pm 2 ); the 
real n k (X) and imaginary kk (X) parts of the complex refrac- 
tive index; and the fraction of the spherical particles C^ sp h- 
The spectral dependence of optical thickness (X)/x k(Xi) 

is also vertically invariant and defined by these parameters, 
while the absolute value of (X) additionally depends on 
the total volume of the aerosol in the atmospheric column: 

c*,v = J2i = i Ni or der to account for vertical 

variability of the aerosol extinction k ex t,& (h, X), the addi- 
tional characteristic Ck (h) was added. The function Ck (h) 
defines the vertical distribution of aerosol concentration and 
the optical thickness of k-th aerosol component in each of 
i-th atmospheric layer is defined as: 

hi 

A r ijk (X) = Tk (X) J c k (h) dh. (16a) 

hi + 1 

The aerosol concentration profile Ck (h) is assumed as a 
Gaussian function normalized to unity, i.e.: 
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ck (h) 


exp - 


(h - h k ,o) 


) hroA 

-/ 

^BOA 


c k (h) dh = 1, (16b) 


where hsoA - Bottom Of the Atmosphere (BOA) height and 
hjOA ~ Top Of the Atmosphere (TOA) height. In the re- 
trieval, the standard deviation characterizing the width of the 
aerosol layer is fixed to <7k = 0.75 km. Therefore, only one 
parameter characterizing aerosol vertical distribution is in- 
cluded into the retrieval: hk , o - the mean altitude of the 
k-th aerosol component layer. The present version of the 
POLDER/PARASOL retrieval code is set to retrieve only one 
aerosol component. The algorithm derives single values of 
complex refractive index and fraction of spherical particles 
for particles of all sizes. Such retrieval assumption is based 
on the earlier AERONET sensitivity studies by Dubovik et 
al. (2000, 2006) that indicated major limitations in discrimi- 
nating between refractive indices and shapes of aerosol par- 
ticles of fine and coarse modes. At the same time, the possi- 
bility of retrieving several aerosol components with different 
complex refractive indices and vertical distributions is also 
assumed (see Table 1) and sensitivity of polarimetric obser- 
vations to multi-component aerosols is planned to be verified 
in future studes. 

In addition, in order to harmonize the radiative trans- 
fer code with the structure and needs of general inver- 
sion approach, several modifications have been implemented. 
Specifically, the modifications were aimed to increase the 
speed of calculations by allowing admissible decrease of the 
accuracy of modeling. The three possible tradeoffs permit- 
ting reduction of computation time without any significant 
loss of retrieval accuracy were identified and implemented. 


3.3.1 Adjustable number of terms in the expansion of 
the phase matrix and in the quadrature of 
directional integration 


The accuracy of radiative transfer calculations strongly de- 
pends on the number of terms M used in the expansion of 
the phase matrix into Legendre polynomials and number of 
terms N used in Gaussian quadrature for zenithal integra- 
tion. The values should satisfy the inequality 4 N — 1 > 2 M 
to retain conservation of energy in the successive order of 
scattering integration. The values M and N should be suf- 
ficiently large to provide accurate calculation. However, the 
larger M and N the longer the calculation time. At the same 
time, the high accuracy of the modeling is not always re- 
quired during the retrieval. For example, studies by Dubovik 
and King (2000) showed that when observations of ground- 
based radiometers are inverted, the successful retrieval can 
be achieved using the approximate and quick calculations 
of the first derivatives. Correspondingly, the retrieval time 
can be significantly decreased because the calculations of 
the first derivatives is the most time consuming component 
of Newtonian’s retrieval algorithms. Following this strat- 
egy, the operational retrieval of aerosol from AERONET data 


relies on analytical single scattering approximation for cal- 
culations of derivatives. The possibility of using this trade- 
off in POLDER retrieval algorithm was tested. It has been 
concluded that using single scattering approximation is not 
sufficient for conducting retrieval from satellite observations. 
Nonetheless, it has been found that the Jacobians estimated 
numerically as finite differences on basis of the full radiative 
transfer calculations implemented with significantly reduced 
values of M and N provide fast and accurate retrievals. This 
approach is used in the present algorithm. 

The utilization of linearized radiative transfer code 
(e.g. see Hasekamp and Landgraf, 2005a) that provides all 
derivatives with respect to aerosol and surface properties in a 
single run would be alternative and promising strategy of ac- 
celerating satellite observation inversion. This strategy was 
not used in the present study because the linearization of ra- 
diative transfer code is a rather complex effort that requires 
significant time investments. The possibility of using this ap- 
proach will be considered in future studies. At the same time, 
a retrieval algorithm relying on the numerical calculation of 
the first derivatives is probably more flexible in practical ap- 
plications. Indeed, in the present POLDER/PARASOL algo- 
rithm the set of retrieved aerosol or/and surface parameters 
can be changed with no modifications in the calculations of 
the first derivatives. If the derivatives are calculated analyti- 
cally, achieving such flexibility could be more difficult. 

3.3.2 Truncation of the phase matrix 

The truncation of the phase function is a technique whereby 
the scattering effects from the sharply increasing forward 
peak of the phase function are calculated separately from 
those of the rest of the phase function, which permits the 
accurate but much faster modeling of diffuse radiation. For 
example, the AERONET aerosol operational retrieval by 
Dubovik and King (2000) employs the discrete ordinate ra- 
diative transfer code by Nakajima and Tanaka (1988) that 
uses the efficient procedure of the phase function trunca- 
tion and provides very fast and accurate calculation of down- 
welling diffuse radiation in moderately thick atmospheres. 
The detailed discussion of different methods of implement- 
ing the phase matrix truncation and their comparison are 
given in the recent paper by Rozanov and Lyapustin (2010). 

Following the ideas of Nakajima and Tanaka (1988) the 
truncation of the phase matrix has been implemented in the 
successive order of scattering code as a part of present stud- 
ies. In the developed version of the PARASOL algorithm the 
use of truncation is optional but recommended. The utiliza- 
tion of the phase matrix truncation allows for decreasing the 
number of terms M in the expansion of the truncated phase 
function and N in the Gaussian quadrature for azimuth inte- 
gration. According to the results of the conducted tests, accu- 
rate PARASOL retrievals can be achieved with the following 
recommended values: M = 21 and N = 10 - for calculating 
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the fit to PARASOL observations; M= 15 and N = 7 - for 
calculating Jacobian matrices. 

3.3.3 Flexible angular representation in the phase 
function 

The successive order of scattering radiative transfer code 
uses the phase function values at /V ang angles correspond- 
ing to the points of Gaussian quadrature. These values are 
provided to the radiative transfer code from the Dubovik et 
al. (2006) software generating aerosol single scattering prop- 
erties. The software provides the phase function for the 
set of fixed scattering angles allowing sufficiently accurate 
modeling of angular variability of scattering. For example, 
in AERONET retrieval the phase function is calculated at 
83 Gaussian points. At the same time, the speed of aerosol 
single scattering modeling depends on the number of the 
scattering angles used. In order to enhance the flexibility of 
the retrieval code the possibility of using the phase function 
at nearly arbitrary set of scattering angels has been imple- 
mented, with the values of the phase function at the required 
Gaussian points derived from the software data by conve- 
nient interpolation. The results of the series of the numerical 
tests have shown that POLDER/PARASOL observations can 
be adequately modeled using set of A ang = 35 selected scat- 
tering angles. 

4 Numerical inversion 

In contrast to the majority of existing satellite retrieval algo- 
rithms, this effort is one of the first attempts to develop an 
aerosol satellite retrieval using statistically optimized multi- 
variable fitting. Such strategy does not rely on the pre- 
assumed classes of potential solutions. Instead the solution 
is sought in a continuous space of solutions under statisti- 
cally formulated criteria optimizing the error distribution of 
the retrieved parameters. The implementation of some ele- 
ments of such a strategy was pursued in the earlier develop- 
ments of satellite retrieval algorithms. For example, the sta- 
tistical optimization of the retrieval solutions was used for in- 
version of MISR observations by Martonchik et al. (1998), in 
the retrieval algorithms proposed by Chowdhary et al. (2002, 
2005) and by Waquet et al. (2007, 2009a) for inverting antic- 
ipated GLORY/APS observations and applied to RSP data, 
in the retrieval approach suggested by Hasekamp and Land- 
graf (2005b, 2007) for applying to multiple-viewing-angle 
intensity and polarization measurements and the retrieval al- 
gorithm developed by Govaerts et al. (2010) for processing 
SEVIRI/ MSG observations. 

Detailed descriptions of inversion methods can be found 
in various textbooks (Tikhonov and Arsenin, 1977; Twomey, 
1977; Tarantola, 1987; Press et al., 1992; Rodgers, 2000; 
Doicu, 2010). However, the textbooks generally provide 
the reviews of well-established inversion procedures and do 


not suggest the reader sufficient explanations as to which 
method and why should be chosen for a particular applica- 
tion. The approach used here is focused on clarifying the 
connection between different inversion methods established 
in atmospheric optics and unifying the key ideas of these 
methods into a single inversion procedure. It follows the de- 
velopments by Dubovik and King (2000), Dubovik (2004), 
Dubovik et al. (2008). The methodology has several orig- 
inal (compare to standard inverse methods) features opti- 
mized for remote sensing applications. As shown in de- 
tailed description by Dubovik (2004), the methodology ad- 
dresses such important aspects of inversion optimization 
as accounting for errors in the satellite observations, in- 
version of multi-source data with different levels of ac- 
curacy, accounting for a priori and ancillary information, 
estimating retrieval errors, clarifying potential of employ- 
ing different mathematical inverse operations (e.g. com- 
paring iterative versus matrix inversion methods), accel- 
erating iterative convergence, etc. The concept uses the 
principles of statistical estimation and suggests a general- 
ized multi-term Least Square type formulation that com- 
plementarity unites advantages of a variety of practical in- 
version approaches, such as Phillips -Tikhonov -Twomey con- 
strained inversion (Phillips, 1962; Tikhonov, 1963; Twomey, 
1963), Kalman filter (Kalman, 1960), Newton-Gauss and 
Levenberg-Marquardt iterations, etc. This approach pro- 
vides significant transparency and flexibility in development 
of remote sensing algorithms for deriving such continuous 
characteristics as vertical profiles, size distributions, spec- 
tral dependencies of some parameters, etc. For example, 
compared to the popular “Optimal Estimation” equations 
(Rodgers, 2000), the multi-term Least Square type formu- 
lation allows harmonious utilization of not only a priori esti- 
mate term but instead, or in addition, using a priori terms lim- 
iting derivatives of the solution (see discussion by Dubovik, 
2004 and Dubovik et al., 2008). This methodology has re- 
sulted from the multi-year efforts on developing inversion 
algorithms for retrieving comprehensive aerosol properties 
from AERONET ground-based observations. 

Two alternative scenarios are proposed for inverting satel- 
lite observations: single-pixel retrieval and multiple -pixel re- 
trieval. The single-pixel retrieval is a conventional approach 
when observations of the satellite instrument over each single 
pixel (e.g. in the case of POLDER/PARASOL the pixel size 
is 5.3 km x 6.2 km at nadir) are inverted completely indepen- 
dently. The multiple-pixel retrieval is a newly suggested ap- 
proach when the observations of the satellite instrument over 
a group of pixels are inverted simultaneously and extra a pri- 
ori constraints on the inter-pixel variability of the retrieved 
parameters is applied. It is expected that applying such con- 
straints will help to improve reliability of the retrieval. 
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4.1 Single-pixel retrieval 

This part of the inverse procedure is adopted, with some 
modifications, from the AERONET (Dubovik and King, 
2000) inversion algorithm. There are two key aspects in the 
algorithm organization: the general organization of observa- 
tionfitting and the a priori data representation. 

4.1.1 Single-pixel observation fitting 

Formally, the retrieval algorithm is designed as a multi-term 
LSM that implements statistically optimum fitting of several 
sets of observations and a priori constraints under assump- 
tion of normally distributed uncertainties (detailed discussion 
is given by Dubovik, 2004). It provides the solution of the 
following system of equations: 

/* = /(«) + A / 

0* = (A a)* = Sa + A (A a) . (17) 

a* — a A a* 


The statistically optimized solution of Eq. (17) defined on 
the base of Maximum Likelihood Method corresponds to the 
minimum of the following quadratic form: 

'h (fl p ) = % (fl p ) + 'I'a (fl p ) + 'I'a («P) ( 18 a) 

= l - ((A fP) T (W/) -1 A /? + K a (fl p ) r a « p 

+ K a (« p - a*) T (« p - a*)). 

In generally non-linear case the minimum can be obtained by 
the following iterative procedure: 

fl p +i = fl P _ fp a o p , (18b) 

where A a p is a solution of the p-th so-called, in statistical 
estimation formalism, Normal System 

A p A fl p = V4 («P). (18c) 


Here, /* is a vector of the PARASOL measurements, A/ is 
a vector of measurement uncertainties, a is a vector of un- 
knowns. The second line in Eq. (17) represents the a pri- 
ori smoothness assumptions used to constrain variability of 
size distribution and spectral dependencies of the real and 
imaginary parts of the refractive index as well as the spec- 
tral dependencies of parameters of surface reflectance model. 
The matrix S includes the coefficients for calculating ra- 
th differences (numerical equivalent of the derivatives) of 
dV(rj)l dlnr, n(Xj), k (Xj), p Q (Xj), k (Xj), 6 (Xj); 0* - vec- 
tor of zeros and A (A a) - vector of the uncertainties char- 
acterizing the deviations of the differences from the zeros. 
Formally, this equation states that all these ra-th differences 
are equal to zeros within the uncertainties A (A a). The third 
line in Eq. (17) includes the vector of a priori estimates a'* 
and A a* is the vector of the uncertainties in a priori esti- 
mates. 

In order to account for the non-negative character of 
the observed radiances and retrieved aerosol (dV(rj)l dlnr, 
n (Xj), k (Xj)) and surface reflectance (p 0 (Xj)) parameters, 
the assumption of log-normal error distribution was used. 
The noise distribution is apparently the most appropriate for 
the positively defined values. The log-normal noise distribu- 
tion implies that the logarithms of the observed positively de- 
fined values are normally distributed. Thus, for convenience 
of formulating the statistically optimized solution of Eq. (17) 
the logarithmic transformation is used for both measured fj 
and retrieved at parameters (see detailed discussions on va- 
lidity of applying this transformation in the publications by 
Dubovik and King, 2000, Dubovik, 2004). Correspondingly, 
the errors A/, A (A a), and A « * are assumed normally dis- 
tributed. Table 2 shows the exact definitions of each element 
of the vectors /* and a. In addition, Table 2 shows the vari- 
ability ranges allowed for each retrieved parameter. 


The matrix in the left side of Eq. (18c) is also known as 
Fisher Matrix A p and the right side of Eq. (18c) represents 
the gradient (a p ) of quadratic form (a p ) in vicinity of 
a p : 


A p = Kp W f -‘ K p + ka A + Ka W a -\ 

(18d) 

V (a p ) = 

(1 Se) 


= Kp W 1 A / p + ka a? + Ka (a p - a*), 


where A/ p = / (a p ) — /*, K p - Jacobi matrix of the first 

derivatives ^ . It should be noted that Fisher Ma- 

trix A p can be considered as so-called Hessian matrix of 
second-order partial derivatives of the quadratic form (a p ) 
(e.g. see Bevington, 1969; Tarantola, 1987). Correspond- 
ingly, Eq. (18c) can be also written as follows: 


(VV4/ (« p )) Aa p = V4/ (« p ). 


(1 8f> 


where VV r 4> (a p ) is the matrix with the elements 


{w r * («%,= 


a—aX 


In the present inversion strategy, all equations are ex- 
pressed via W - weighting matrices that defined by dividing 
the corresponding covariance matrix C by its first diagonal 
element e 2 , i.e. W = (1 /e 2 )C. This formal transformation of 
the inversion equations allows rather transparent interpreta- 
tion of Lagrange multipliers y - parameters determining the 
contributions of a priori terms into solution. Following the 
concept proposed by Dubovik and King (2000) the Lagrange 
multipliers y a and y& are defined as: 


Ka = s}/s\ and Ka = ej/s a , (19a) 

where ej , s\ and are the first diagonal elements of corre- 

sponding covariance matrices Cf, Ca and C a . Note, that if 
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Table 2. Definition of the measurement vector /* and the vector of unknowns a. 


/* - vector of measurements: 

{/* } . = ln(/ (@y ; A,;)), where I (@y ; A,;) is total radiance observed by POLDER/PARASOL 

{/*}.= ln(P (©y; A./)), where P (0^; A,/) = ’ is degree of linear polarization 

obtained from POLDER/PARASOL observations 


a - vector of unknowns 


Notation 

Definition 


Variability limits 

AEROSOL parameters: 

Uy 

(«vh= 

..., N r 

0.00001 < 

ay c 

ay c = ln(C v ) 


0.001<C V < 5.0 (pm 3 /pm 2 ) 

a sph 

a sph — l n (L S ph) 


0.001<C sph < 1.0 

a n 

{a n )i = ln(w(A,-)), i = 1, .. 


1.33 <n(Xi) < 1.6 

a k 



0.0005<fc(A,i)<0.1 

Surface BRDF and BPDF parameters: 

fl brdf,l 

{«brdf,l }/ =ln(po(^/», i : 

= h...,N x 

0.001 <po(ki)< 0.7 

«brdf,2 

{«brdf, 2 };=ln(K(A;)), i = 

-- L N x 

0.1 <k(X() <1.0 

«brdf,3 

{^brdf, 3 } / — l n (l+$(7/)). 


—0.5 < 0{Xi) < 0.5 

tt brdf,4 

{ fl bidf,4}/ = ln(/* 0 (k*)), i 


0.001</z 0 (A-/) < 0.7 

tt bpdf 

{a hvd f} i =\n(B(X i )),i = 

h...,N x 

0.01<#(A;) < 10.0 


one can assume a simplified situation when the POLDER ob- 
servations /*, a priori estimates of differences (A a)* and a 
priori estimates of parameters a * are independent and have 
the same accuracy, i.e. C f = ejl( NfxN{ ), C a =s 2 a I( N axNa) 
and C a = sll(N a xN a ), then corresponding weighting matri- 
ces are simply equal to the unity matrices: Wf = I(# f xiVf)» 
W A =I( 7 V a xAa) an d Wa=I(Af a xA a )- 

In addition, as discussed by Dubovik (2004), straightfor- 
ward increasing of Nf by adding redundant very similar ob- 
servations is not necessarily beneficial for the retrieval. For 
example, from a practical viewpoint, adding many observa- 
tions with identical or nearly identical observation conditions 
does not bring new information and does not improve the ac- 
curacy of the solution. However, formally even such addition 
of observations artificially increases the impact of the mea- 
surement term on the solution compare to a priori second and 
third terms in Eq. (18). Therefore, the Lagrange multipliers 
y are scaled as: 


YA 


Nfsj 

Nas 2 a 


and y a = 


Nfsj 

N,s 2 ’ 


(19b) 


As suggested by Dubovik (2004) such scaling reflects the in- 
evitable decrease of accuracy of single observation s 2 if the 
number of observation Nf increases excessively. The above 
equation is written under an assumption that increasing the 
number of measurements in the coordinated set of remote 


sensing observations inevitably results in a decrease of the 
accuracy of each single measurement in this observation set. 
For example, if a satellite sensor takes one single observa- 
tion, the expected variance of measurement error is s 2 N . If 
the same sensor makes Nf space- and/or time-coordinated 
observations the variance of the error in each single obser- 
vation increases by the factor Nf, i.e. e 2 N ~ NfS 2 v This in- 
crease can be explained by the fact that the consistency of 
the Nf coordinated observations should be assured by con- 
trolling Nf relations between the Nf observations. The con- 
trol of each of those relationships introduces a random error 
s 2 n , correspondingly the error variance of a single measure- 
ment in Nf dimensional observation increases in Nf times. 

The coefficient t v < 1 in Eq. (18b) is adjusted to provide 
monotonic decrease of 4* (« p ), i.e. 

V (a p+1 ) < 'I' (« p ). (20a) 

If all assumptions are correct, the minimum value of the 
above quadratic form can be theoretically estimated as fol- 
lows: 

* (a) (Nf + N a + JVa* - JVa) Sf . (20b) 

Note that the minimum value of W (a p ) relates to e 2 because 
the weighting matrices were used instead of using directly 
the covariance matrices. Once the value of measurement er- 
ror is known s 2 , Eq. (20b) can be used to verify the consis- 
tency of the retrieval. Specifically, the inability to achieve 
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the above minimum can indicate the presence of unidentified 
biases or inadequacy in the assumptions made. 

It should be noted that the control of “measurement resid- 
ual” fi'f (a p ) (the first term of quadratic form in Eq. (18a)) 
is a very useful tool for diagnostics of the retrieval dynam- 
ics. Specifically, the final value of (a?) should be close 
to the level of the expected measurement noise. Indeed, if 
the algorithm has found the right solution the value of the 
total residual (a p ) should be rather small and determined 
mainly by the random errors of observations. The contri- 
butions of the a priori residual terms in Eq. (18a) should 
not be significant, because generally the weights of a pri- 
ori terms 4^ (« p ) and 4> a (a p ) are minor compared to the 
weight of the “measurement residual” term (a p ). How- 
ever, at early iterations when the solution approximation is 
very far from the solution, 4q (a p ) is dominated by lineariza- 
tion errors and has the value much higher than the level of 
the expected measurement noise. Therefore, since accuracy 
of a priori data is independent of the iteration, the weight of 
the a priori term should be increased. Correspondingly, this 
additional enhancement of the a priori data impact on the so- 
lution improves the convergence of non-linear fitting. For 
example, in developed POLDER/PARASOL algorithm, fol- 
lowing Dubovik and King (2000) and Dubovik (2004), the 
strength of a priori constraints is adjusted dynamically as a 
function of the measurement residual 4>f (a): 


K) 


(« p ) 

(N f - Ay ’ 


(20c) 


where 4q (« p ) reflects the accuracy of the POLDER observa- 
tion fit at the p-th iteration and provides an indication of how 
the iterations converge to the solution. For example, at the 
last iteration, when the solution estimate « Plast is expected to 
be in a small vicinity of the actual solution, the value of the 
residual 4>f (a) of the measurement fit can be estimated as: 


q/ f ( fl Pb8t) « (N f - A a ) £f. 


(20d) 


Using sj (« p ) Eq. (20c) adjusts the values of the Lagrange 
multipliers ya and y a in Eq. (18) and enhances the contri- 
bution of a priori constraints at the earlier iterations. As 
suggested by Dubovik (2004), this dynamic determination 
of a priori constraints improves convergence of non-linear it- 
erations analogously to Levenberg-Marquardt formulations. 
At the same time, in contrast to the original Levenberg- 
Marquardt method, the idea of enhancing constraints on 
the solutions at earlier iterations in the formulations by 
Dubovik (2004) is included harmoniously within the frame- 
work of united statistical estimation approach. For exam- 
ple, if no smoothness constraints are used (i.e. the smooth- 
ness terms in the right sides of Eqs. (17) and (18) are elimi- 
nated) and no a priori estimates a * are available, then one can 
assume a*=a p , Eq. (18) become equivalent to Levenberg- 
Marquardt formulations (see details in Dubovik, 2004). 


Thus, the inversion procedure described above is driven 
by the limited set of the input characteristics that includes 
the weighting matrices W ... and corresponding variances e . 
The vector of POLDER/PARASOL measurements includes 
2 components: 



(21a) 


where index “I” denotes total reflectance observations and 
index “P” denotes the observations of the degree of linear po- 
larization (see Table 2). Since the total radiance and degree 
of linear polarization are measured with different accuracy, 
the covariance matrix assumed for logarithms of measure- 
ments /* has the array structure: 



(21b) 


Correspondingly, the weighting matrix Wf is defined as: 



(21c) 


where I... denote the unity matrices of corresponding dimen- 
sion and yp is a ratio of the total and polarized reflectance 
variances: 


Yp = fip/fii 2 - (2 Id) 

The variance of the errors in measurements of total re- 
flectance is expected at the level of 2 % relative to the mag- 
nitude of observed radiance /, i.e. s\ = A (In/) & =0.02. 

The variance of the errors in measurements of the degree of 
linear polarization is expected at the absolute level of 1 %, 
i.e. AP = 0.01 and s p = A (InP) ss = qji. Equation (21b) 
assumes the simplest structure of the covariance matrix when 
intensity and polarization observations are equally accurate 
for all spectral channels and angles of observations. If more 
detailed information about covariance matrix is available it 
can be trivially integrated into Eq. (21). 

In contrast to the numerous remote sensing algorithm re- 
lying on a priori estimates as suggested by Rodgers (2000), 
the a priori estimates a * were not used at all, i.e. oo 
(although in practical application of the algorthm to satellite 
observation the using of, at least for land surface parameters, 
the climatological values with adequately estimated vari- 
ances e\ can be beneficial). Thus, in contrast to the method- 
ology by Rodgers (2000), the retrieval was constrained here 
using exclusively the a priori smoothness constraints as dis- 
cussed in the follow on Section. 


4.1.2 A priori smoothness constraints in single pixel 
fitting 

The vector a includes several components: 
a = (#v & n #k *7sph a Vc #brdf, 1 #brdf,2 #brdf,3 #bpdf) 5 (22) 


www.atmos-meas-tech.net/4/975/2011/ 


Atmos. Meas. Tech., 4, 975-1018, 2011 



992 


O. Dubovik et al.: Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties 


where a w , a n , dk, <2 S ph denote the components of the vector a 
corresponding to the dV(rj)/dlnr, n (kj), k (A/) and C, <2his 
equal to the logarithm of mean altitude of the aerosol layer 
h a . The element ciy c is the logarithm of total volume concen- 
tration, while a w includes logarithms of values dV(r)/d\nr 
normalized by total volume concentration. This modifica- 
tion (compare to AERONET algorithm) allows decoupling 
the retrieval of the total amount of aerosol from the shape 
of the size distribution. This decoupling of the size distri- 
bution parameter is essential for satellite imager algorithm, 
since it allows applying different constraints on the shape of 
size distribution and aerosol loading. The three components 
«brdf,i> «brdf,2 and «brdf,3 that are related to the logarithms 
of spectrally dependent parameters p Q (k), k (k) and 0 (k) 
of employed RPV model. Note, that here the “hot spot” pa- 
rameter ho (k) that is present in Eq. (12) is not included in the 
vector d, because it is not retrieved near-backscattering direc- 
tion nor observed (it is significant only in a narrow range of 
scattering angles around the backscattering direction). At the 
same time, ho (k) is retrieved if the “hot spot” is observed. 
In such situation it is included in the retrieval similarly to 
other BRDF parameters (p Q (k), k (k) and 0 (k)). The vec- 
tor « bpdf includes the logarithms of the spectrally dependent 
free parameter B ( k ) of the employed Maignan et al. (2009) 
model. The detail description of each element of the vector 
d is given in Table 2. 

The a priori smoothness constraints are applied in the 
POLDER/PARASOL algorithm on several different compo- 
nents of the vector d differently. For example, for d given by 
Eq. (22), the matrix S has the following array structure: 

/S v 00000 0 0 0 0 

OSJOOO 0 0 0 0 

0 0 S k 0 00 0 0 0 0 

0000000 0 0 0 

0 0 0 000 0 0 0 0 

0000000 0 0 0 

0 0 0 OOOSbrdfi 0 0 0 

0 0 0 0 0 0 0 Sbrdf,2 o 0 

000000 0 0 S brd f,3 0 

^ 0 0 0 000 0 0 0 Sbpdf 

where the corresponding matrices S... have different dimen- 
sions and represent differences of different order (3 for size 
distribution, 1 for n (k), 2 for k (k), 2 for p Q (k) 1 for k (k), 
0 (k) and B ( k ). The lines in Eq. (23) corresponding to a h, 
ay c and <2 sp h contain only zeros because no smoothness con- 
straint can be applied. It should be noted that for making 
formulations more transparent only one coefficient ya was 
shown in Eq. (18). However, the actual algorithm uses 7 (or 
even 8 if h o (kj) is retrieved) multipliers Yaj- The errors 
A (Ad) are assumed independent for different components 
of the vector (A d)* and the smoothness matrix in Eq. (18) 
has the following array structure: 



y A ,i £2i 0 0 000 0 0 0 0 

o KA,2^2 o 000 0 0 0 0 

0 0 y A , 3 fl 3 000 0 0 0 0 

0 0 00000 0 0 0 
000 000 0000 
000 000 0000 
0 0 0 0 0 0 y A> 4 ^4 o 0 0 

0 0 0 0 0 0 0 Ya,5 ^5 0 0 

0 0 0 0 0 0 0 0 y A ,6 ^6 o 

0 0 0 0 0 0 0 0 o y A> 7 ^7 

where ft/ = Sj Wr 1 S/ uses the derivative matrices S; 0 = 1, 
...,7) S v , S n , Sk, Sbrdf, i , Sbdrf.25 Sb r( jf, 3 , S^p^f. If it is assumed 
that the covariance matrices of errors A ( Ad) have the struc- 
ture Ca,/ = £a i l (N A /xVa i) f° r each components of (Ad)*, 
i.e. weighting matrices are Wa,/ = I(/v a i*N A ,•)• Correspond- 
ingly, the quadratic form 4 >a 0* p ) in Eq. (18) can be written 
as the following sum: 

2 4< a («p) = )/ A («P) r S2aP= (24) 

= E 2 ' l ’ A - i ( flP ) = E XA.i {a]) T S2, a P. 

i= 1,...,7 

where y A ,/ =e f / e A,r 

The utilization of the smoothness constraints for a sin- 
gle retrieved function y(xi) was originated in the papers by 
Phillips (1962), Tikhonov (1963) and Twomey (1963). Al- 
though application of the smoothness constraints is usually 
considered to be an implicit constraint on derivatives, in these 
original papers and most of follow-on studies the solution d 
was constrained by minimizing ra-th differences A m of the 
vector d components: 

A 1 = dj + 1 — dj, (m = 1), (25a) 

A 2 = a i+ 2 - 2 flj+i + dj (m = 2), 

A 3 = a i+ 3 - 3 d i+ 2 + 3 d i+ i - dj, (m = 3). 

The corresponding “smoothness” matrix ft = (Sk) r (S^) was 
defined using S( m ) matrices of ra-th differences S( m ) 
(i.e. A m =S ( m )d). For example, S( 2 ) (ra =2) is: 

/ 1 -2 1 0 ... \ 

0 1-2 1 0 ... 

S 2 = 0 0 1 -2 1 0 ... . (25b) 

v 0 1-21 j 

The present development follows the concept of Dubovik 
and King (2000) and Dubovik (2004) that considers smooth- 
ness constraints explicitly as a priori estimates of the deriva- 
tives of the retrieved characteri sticy (xj ) . The values of ra-th 
derivatives g m of the function y(x ) characterize the degree of 
its non-linearity and, therefore, can be used as a measure of 
y(x) smoothness. For example, smooth functions y(x), such 
as, a constant, straight line, parabola, etc. can be identified 
by the ra-th derivatives as follows: 
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gi CO = dy (x)/dx = 0 => yi (x) = C; (26) 

g 2 (x) = d 2 y ( x)/dx 2 = 0 =>> y 2 (x) = B x + C; 

g 3 (jc) = d 3 y (x)/d* 3 = 0 =)> y 3 (x) = A x 2 + B x + C. 

These derivatives can be approximated by differences be- 
tween values of the function a; = y(x/ ) in A a discrete points 
Xi as: 

dy (**') _ A 1 y (x;) y (x t - + A x,-) - y(x/) 

« = (27a) 

ax Ai Xi Ai Xj 

_ y (x i+ 1 ) - y(xj) 

Ai Xi 

d 2 y (xj») _ A 2 y (x t -) 
r/x 2 A 2 (x;) 

_ A l y(x i+ i)/Ai(x i+ i) - A 1 y (x,-)/Ai (x/) _ 

(Al Xi + Al X; + i)/2 

y(x t -///) _ A 3 y(xj) 
dx 3 A 3 (x/) 

_ A 2 y (X; + i)/A 2 (x/ + i) - A 2 y (x t -)/A 2 (x,-) _ 

(A 2 (x,0 + A 2 (xi+i))/2 

where 

Ai (x/) = Xi+i-x,-; 

A 2 (x/) = (Ai (x/) + Ai (x,-+i))/2; 

A 3 (x/) = (A 2 (x,-) + Ai (x/+i))/2; 
x,-/ = Xi + Ai (x/)/ 2 ; 


Xi" = X| + (Ai(X|-) + A 2 (x/))/2; 

x/'" = x/ + (Ai(x/) + A 2 (x/) + A 3 (x/))/2. 


The corresponding matrix of m-th derivatives S gj ( m ) 
(i.e. = Sg,( m ) a) can easily be defined using Eq. (25). For 

example, S g> ( 2 ) ( m = 2 ) is: 


S «.(2) = 


(27b) 


/ 2 -2 2 

Ai (Ai + A2) (Ai A2) A2 (Ai + A2) 

0 2 = 2 

^ 2 « (a 2 -a 3 ) 

0 0 ? 

u u A 3 (A 3 + A 4 ) 


0 

2 

-2 

(A3 a 4 ) 


0 


...\ 


2 

A 4 (A3 + A 4 ) 


0 


V ) 

where A,-=Ai (x,-) =x,-+i — x,-. One can see that S g> ( 2 ) 
has significantly more complex structure than the conven- 
tional definition of smoothing matrix by Eq. (25). At the 
same time for the special case of uniformly spaced ordinates 
Ai = A = const the matrices S« ( OT ) and S( m ) have straightfor- 
ward relation: S gj ( OT ) = A -( 2 m) S( m ). 


In difference with Eq. (25), Eq. (26) allows for apply- 
ing smoothness constraints in more general situations when 
Ai (x;) / const. For example, in present algorithm there is 
a number of spectral parameters that are functions of X and 
the algorithm deals with their values defined for each spec- 
tral channel A/ . Obviously the (k ; + 1 — A./) / const and using 
standard definition of differences by Eq. (25) for smoothing 
spectral parameters (e.g. n(X), p Q (X) etc.) is not completely 
correct. Applying the limitations on the derivatives defined 
by Eq. (27) is more rigorous if no significant changes of 
derivatives of y(x) are expected for different ordinates x/. 
Although using Eq. (27) leads to a loose of transparency in 
definitions of matrices S( m ), generating those matrices on al- 
gorithmic level is rather straightforward and defining of La- 
grange parameters ya is more logical. Therefore, in present 
algorithm, the smoothness constraints are applied to limit di- 
rectly the numerical equivalents of m-th derivatives given by 
Eq. (25), by a priori assumption that m-th derivatives are 
equal to zeros with some uncertainty (g* = 0 * = 0 + A g ). 

In order to estimate correctly the strength of a priori con- 
straints Dubovik and King (2000) and Dubovik (2004) sug- 
gested to relate the strength of a priori smoothness con- 
straints (i.e. the uncertainty s\ • in a priori known deviations 
of derivatives from zeros in second line of Eq. (17) to the 
integral norm of the derivatives defined as: 



where ^characterizes the smoothness of the physical contin- 
uous function y;(x): the closer b t to zero the more smooth 
the function y; (x). Indeed, there is a direct relation of this 
integral to the quadratic forms T'a,; (a v ) in Eq. (24) to the 
integral norm of the derivatives given by Eq. (28). Namely, if 
the values of function yt (x) are known at A/ discrete points 
Xi, one can estimate the norm of the m-th derivatives: 



rr r 1 

where Q g j = (S g j) W- S g j, S g j are the matrices of the 
coefficient for estimating derivatives g m j (x 7 ) as shown 
in Eq. (27), W,- is diagonal weighting matrix with the el- 
ements {Wi}jj = l/(A mA Xj)). Therefore, the multiplier 
YA,i =£f /e\ j in the minimized quadratic form of Eq. (24) 
can be directly related to the known values of bj as follows: 


YA,i 



(30) 


where 


m 


max 



d m y? s (x) 
d™x 


dx. 


(31) 
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Fig. 5. The general structure of the numerical inversion data flow implemented for the retrieval of aerosol over single pixel. 


Here y* s (x) denotes most unsmooth real function yt(x). 
Based on this definition, the smoothness a priori condition 
imposed by the second line in Eq. (17) allows any function 
yi(x ) to be solution that has norm of ra-th derivatives b™ 
smaller than i n the framework of used here sta- 

tistical approach this definition means the second system in 
Eq. (17) defines the derivatives equal zeros: g* = 0* = 0 + 
A i tg . The errors A g are assumed independent for each z-th 
(i = 1, 2, ..., 7) component of vector a (see Eq. 22) have the 
following diagonal covariance matrix with the elements of 
the diagonal defined as: [C g j} .=s 2 g ■ (A m>l - (xj)). This as- 
sumes that the deviation of estimates g m j (v ; ) from zero has 
variance s 2 •/( A m ,z (*/))• This variance is chosen inversely 
proportional to the coefficient A m j (xj) calculated as shown 
by Eq. (27). This inverse proportionality to A m j (xj) re- 
flects the fact that the deviations of g m ,i(xj ) from zeros can 
be stronger at smaller scales. This is essential in general case 
when Ai (xj) / const. Thus, following this consideration we 
define the strength of a priori constraint applied on each re- 
trieved function (dV(r)/d\nr, n (X), k (X), p Q ( X ), k (X), 0 (X) 
and B ( X ) by calculating values of (/?"*) max using most un- 
smooth (variable) examples of corresponding physical func- 
tions, i.e. most unsmooth size distributions, spectral depen- 
dencies of refractive indices and BRDF parameters. Table 3 


shows the orders of the derivatives used for smoothness con- 
straints and the corresponding values of Lagrange multipli- 
ers YA,i- The general structure of the numerical inversion 
data flow implemented for the retrieval of aerosol over single 
pixel is summarized in Fig. 5. 

4.2 Multiple-pixel retrieval 

In contrast with most satellite retrievals, the algorithm de- 
veloped here does not implement the measurement fitting for 
each single pixel independently. Instead, the fitting is imple- 
mented for a group of pixels and is constrained by the extra a 
priori limitations on inter-pixel variability of aerosol and/or 
surface reflectance properties. This approach improves the 
stability of satellite data inversions because the information 
content of the reflected radiation from single pixels is some- 
times insufficient for a unique retrieval of all retrieved param- 
eters. For example, deriving aerosol properties over bright 
land is known as an extremely difficult task. On the other 
hand, if the surface reflectance and aerosol properties happen 
to be the same for a certain time period or over some area, 
one can trivially achieve higher redundancy in the retrieval 
by applying single pixel algorithm to the observations col- 
lected from such multiple pixels. This fundamental tendency 
can be formulated as a priori known statistical limitations that 
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Table 3. The finite differences used for smoothness constrains and 
the correspondent values of Lagrange multipliers YA,i- 



Order of finite 

Values of 


differences 

YA,i 

AEROSOL: 

dV(rj)/d\nr (i = 1, ..., N T ) 

3 

0.005 

n(Xi)(i = l,...,N x ) 

1 

0.1 

k(Xi)(i = l,...,N x ) 

2 

0.01 

Surface BRDF: 

Po(Xi) (i = 1, ..., N x ) 

3 

0.005 

K(Xi) (i = 1, ..., Nx) 

1 

0.1 

0(Xi)(i = l,...,Nx) 

1 

0.1 

ho(Xi) (i = 1, ..., Nx) 

2 

0.01 

Surface BPDF: 

B(Xi)(i = l,...,Nx) 

1 

0.1 


can serve as an extra constraint making multi-pixel retrieval 
more robust and reliable. Dubovik et al. (2008) discussed us- 
ing this concept in a frame of statistical estimation formalism 
for improving global aerosol inverse modeling retrievals. In 
principle, this approach allows for applying smoothness on 
4-D (temporal and spatial: x-y-z) variability of the retrieved 
characteristics, while for inverting POLDER/PARASOL im- 
ager observations only 3-D constraints (temporal and hor- 
izontal: x — y) are used. Methodologically quite similar 
concept was used by Quaife and Lewis (2010) for imposing 
temporal (1-D) smoothness constraints on surface reflectance 
variability for inverting linear kernel-driven BRDF models 
from MODIS observations. It should be noted that applica- 
tion of the multi-pixel approach is discussed here only as a 
tool for applying extra spatial and temporal constraints on the 
retrieved aerosol and surface parameters. At the same time, 
this approach when a group of pixels is inverted together is 
also convenient for rigorous accounting for cross-pixel cou- 
pling in this forward modeling of reflected radiances. This 
coupling appears due to the multiple light scattering effects 
and it introduces some dependence of reflected radiances 
over observed pixel on the properties of surface and aerosol 
over neighboring pixels. In satellite remote sensing this is 
known as adjacency effect that may introduce some addi- 
tional errors in the aerosol retrieval. This effect is significant 
for high-resolution observations and becomes negligible for 
satellite images with resolution ~2km and larger (e.g., see 
Tanre et al., 1981; Kaufman, 1982; Lyapustin and Kaufman, 
2001). Since POLDER/PARASOL observations have single 
pixel resolution of ~5.3kmx 6.2km, the adjacency effect 
was not accounted for in our algorithm. In addition, accurate 
accounting for this effect requires complex and time consum- 
ing 3-D radiative calculations that are difficult to adapt to the 
requirement of the operational satellite retrieval algorithms. 


4.2.1 Multiple-pixel observation fitting 

In order to make the equations more compact Eq. (17) de- 
fined for i- th single pixel can be denoted as follows: 

f /* = fi (a) + A f t 

0* = Sf at + A (A ai ) => /* = f t («/) + A /,-. (32a) 
[ a* = at + A a* 

Then, the inversion of multi-pixel observations can be con- 
sidered as a solution to the following combined system of 
equations: 

ft = fi (fli) + A f x 
ft = f 2 ( a 2) + A f 2 

ft = f 3 (<*3) + A f 3 

(32b) 

0* = S x a + A (A* a) 

0* = S y a + A [A y a) 

0* = S t a + A (A t a) 

where index (i = 1 , 2, 3, ...) denotes the number of each 
single pixel. Correspondingly, the total vector of unknowns 
a is composed from the vectors of unknowns of each i- th 

pixel: 



The matrices S x , S y and S* include the coefficients for 
calculating ra-th differences (numerical equivalent of m-th 
derivatives) of spatial ( x—/y — ) or temporal (t—) inter-pixel 
variability for each retrieved parameter a k characterizing 
dV(rj)/d\nr, n (Xj), k (Xj), Po (Xj), k (Xj), 6 (Xj); 0*, 0J, 
0* - vectors of zeros and A (A x a), A (A y a) and A (A t a) 
- vectors of the uncertainties characterizing the deviations of 
the differences from the zeros. 

The solution of the multi-pixel system can be implemented 
using formally the same sequence of the operations as shown 
in the single pixel case. However, the minimized quadratic 
form T' (« p ), its gradient V*P (« p ) and Fisher Matrix A p 
formulated for inversion of combined multi-pixel Eq. (32b) 
would be defined via corresponding single-pixel terms: 

^pixels \ 1 

y; W; (a p ) I + - (a ,r ) r Winter « p , (34a) 

Alp 0 ... 0 \ 

0 A 2 , p ... 0 

0 0 ... a n p J 

■ V 4/1 (« p ) 

V4/(« p )= Vx M« 2 ) + Si inter a p , (34c) 
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Table 4. The types of the finite differences and the correspondent values of Lagrange multipliers y/±j used in the algorithm for applying 
inter-pixel smoothness constrains. 



Order of finite differences 

Values of y^y 


horizontal (x, y) 

temporal (t) 

horizontal (x, y) 

temporal ( t ) 


continuity 

continuity 

continuity 

continuity 

AEROSOL: 

C v 

2 

- 

0.01 

- 

dV(ri)/d\nr (i = 1, ..., N r ) 

1 

- 

0.1 

- 

6"sph 

1 

- 

0.1 

- 

K 

1 

- 

0.1 

- 

n(Xi)(i = l,...,N x ) 

1 

- 

0.1 

- 

k(Xi)(i = l,...,N x ) 

1 

- 

0.1 

- 

Surface BRDF: 

Po(A f ) (i = 1 , TVa.) 

- 

3 

- 

0.1 


- 

1 

- 

0.1 

0(ki)(i = l,...,N k ) 

- 

1 

- 

0.1 

h 0 (X i )(i = l,...,N k ) 

- 

2 

- 

0.1 

Surface BPDF: 


- 

1 

- 

0.1 


where Tfi (« p ), VTfi (« p ) and A/ ;P are defined for N single 
pixels. The smoothing inter-pixel matrix filter is defined as: 

Sinter = Yx Sj S x + y y S T y S y + y, Sf S t . (34d) 

It is easy to observe that if inter-pixel matrix Winter is elimi- 
nated the solution of the multi-pixel system of N pixels can 
be considered as solution of N independent single pixel sys- 
tems. However, the matrix Winter has non- zero non-diagonal 
elements and VT> (« p ) is not equal to a simple sum of 
VTfi ( a ?). Therefore, the solution of the multi-pixel system 
of N pixels is not equivalent to the solution of N independent 
single pixel systems. At the same time, as shown in the Ap- 
pendix A the inter-pixel smoothness matrix Winter is sparse 
and contains many zeros. This fact allows multi-pixel con- 
straints without a significant increase of computation time. 

Also, the values of Lagrange parameters y x , y y and y t are 
defined using similar concepts to those discussed for apply- 
ing smoothness constraints on variability of aerosol and sur- 
face parameters in single pixel inversion scenario as shown 
in Eq. (29). The only difference is that in Eq. (29) the func- 
tion yi (x) denotes the variability of each single parameter of 
aerosol or surface over x-, y- or t -coordinates. Note, that 
while Eq. (34) suggests the same strength of all constraints 
(that is done for clarity of formulation), the variability limi- 
tations can be very different for each single parameter. This 
is achieved by rather straightforward re-scaling of smoothing 
matrices corresponding to each single limitation (i.e. smooth- 
ing variability of each at over each of x-, y- or t -coordinates 
can be different). Also differences/derivatives of different 


order can be used in each single constraint. This strategy is 
fully implemented in the current algorithm. 

Table 4 outlines the details of the inter-pixel constraints 
employed in the developed algorithm. Figure 6 shows the 
design of multi-pixel inversion algorithm as a rather straight- 
forward generalization of single-pixel inversion. 

The diagram in Fig. 7 emphasizes the similarity of multi- 
pixel inversion with A p i xe i single-pixel retrievals. The dif- 
ference is only in the additions of the inter-pixel smoothing 
terms into the /V p j xe i pixels joint Normal System. The matrix 
Winter defining the inter-pixel smoothing terms is very sparse 
and transparent (see Appendix A), therefore computer time 
required for definition of these smoothing terms is insignif- 
icant compare to other calculations performed in the inver- 
sion. Another difference is that the multi-pixel approach re- 
quires solving the A p j xc i pixels joint Normal System that can 
have rather high dimension. At the same time, the Fisher 
matrix A p of this system is also very sparse and a number of 
numerical tools are available for optimizing the solution of 
linear system with such structure. 

4.2.2 Inter-connection between neighboring groups in 
multiple-pixel fitting 

The previous Section described the simultaneous inversion 
of a group of pixels. The proposed multi -pixel strat- 
egy uses expected continuity in pixel-to-pixel variability 
of aerosol and land surface properties as an additional re- 
trieval constraint. At the same time pixel group has limited 
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The concept of multi-pixel retrieval 



X-Variability Constraints 


Fig. 6. The diagram illustrates the approach for retrieval aerosol 
from a multiple pixels of POLDER/PARASOL. 


size: N x x N y x N t . One can naturally expect the conti- 
nuity/similarity of aerosol and surface properties observed 
the edges of this pixel group and with the properties in the 
pixels-neighbors located just outside of the N x x N y x N t 
pixel group. These additional constraints enforcing continu- 
ity of retrieved properties between different inverted groups 
of pixels can be rather logically added to the inversion for- 
mulations described in the preceding Sections. As shown 
by detailed derivation in Appendix B, the additional condi- 
tions of continuity of aerosol and surface properties with the 
values of corresponding parameter in the neighborhood of 
the inverted pixel group can be implemented by adding ex- 
tra terms in Eq. (34). Specifically, the minimized quadratic 
form ^ (« p ), its gradient VT* (a p ) and Fisher Matrix A p 
will include the extra-terms fledge (« p ), VT'edge (« p ) and 
(Ap)edge accordingly. As demonstrated in Appendix B, these 
additional terms are rather transparent and can be calculated 
with very minor computational effort. At the same time, this 
way of equation organization provides rather powerful addi- 
tional tools constraining the retrieval. Indeed, the inversion 
of satellite observations over large geographical areas and ex- 
tended time periods cannot be implemented simultaneously 
due to natural limitations of computer resources. Instead, the 
large records of satellite observations can be inverted sequen- 
tially by small pixel parcels/groups of limited N x x N y x N t 
size. Adding intergroup constraints makes such sequential 
retrieval nearly equivalent to simultaneous inversion of all 
data. The only difference is that the retrieval for each small 
pixel group would not benefit from the information contained 
in the observations over the pixels inverted in subsequent re- 
trieval acts. At the same time, the retrieval would fully bene- 
fit from all preceding retrievals. Thus, taking into account 
that spatial and temporal correlations of both aerosol and 
surface properties have rather limited range, one can expect 


that sequential pixel parcel-by-parcel retrieval with the inter- 
parcel continuity constraints should produce the extended 
fields of the retrieved parameters. The high consistency of 
those fields are enforced by unified inter-pixel constraints. 

It is worthwhile mentioning the analogy of the sequen- 
tial parcel-by-parcel retrieval with the Kalman Filter type 
retrieval sequences. For example, if the “edge” inter-pixel 
constraints are implemented using only first differences, and 
if they were applied only in the time dimension (when 
the retrieved values of unknowns in the given time are 
used to constrain the solution in consequent time moment), 
then the technique described above in this Section becomes 
full equivalent of known Kalman Filter retrieval technique 
(Kalman, 1960). 

4.2.3 Assuming common parameters for several 
different pixels in the retrieval 

The multi-pixel retrieval scheme suggested in the above Sec- 
tions takes advantage of the limited spatial variability of 
aerosol or land surface reflectance for different pixels or tem- 
porary for the same pixel. However, there are some situa- 
tions when it is reasonable to assume that some parameters 
do not change at all from pixel to pixel. For example, in 
the retrieval algorithms developed for geostationary observa- 
tions, the land surface reflectance parameters and most of the 
aerosol parameters are assumed diurnally constant (e.g. Gov- 
aerts et al., 2010). Moreover, even in the retrievals of aerosol 
from polar-orbiting observations (e.g. POFDER, MODIS), 
neglecting surface property variations on small time scales 
can be a reasonable assumption. Therefore, a possibility of 
assuming the inter-pixel invariant parameters can be a useful 
option for constraining the retrieval and verifying the impor- 
tance of neglected variations. 

Inclusion of such assumptions, in principle, can be imple- 
mented by a significant increase of the corresponding Fa- 
grange multiplier in Eqs. (34) and (B17) (in Appendix B). 
This would enforce very high correlation between corre- 
sponding retrieved parameters and make them practically 
equal. However, such enforcement of the full dependence of 
the retrieved parameters has serious deficiencies. Indeed, the 
increase of y t would result in a strong increase in the smooth- 
ness matrices Winter and £2 e dge into the total Fisher matrix 
in Eqs. (34) and (B18) (in Appendix B). Correspondently, 
making contributions of Winter and £2 e dge dominant would 
doubtlessly produce degenerated Normal systems. There- 
fore, if one can expect that properties of aerosol or surface 
in different pixels do not change, the retrieval of a single 
group of those constant parameters for many pixels seems 
more logical than the retrieval of several groups of param- 
eters (one group per each pixel) forced by the enhanced 
smoothness constraints to have very close values in differ- 
ent pixels. For example, this retrieval approach leads to a 
significantly smaller number of simultaneously retrieved pa- 
rameters that essentially simplifies the numerical inversion. 
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Fig. 7. The diagram illustrating the numerical inversion data flow implemented for the simultaneous retrieval of aerosol properties over 
multiple pixels. 


However, if the algorithm is based on this straightforward 
strategy it looses flexibility, because practical implementa- 
tion of the algorithm allowing easy increase or decrease in 
the number of fixed parameters is logistically very challeng- 
ing. Therefore, here we suggest another approach. The ap- 
proach uses the general formulation of the inversion strategy 
via inter-pixel smoothness constraints even for retrieval sce- 
narios when some derived parameters are fixed to be the same 
for different observations (e.g. assuming that some parame- 
ters of BRDF are intra-day independent over the same pixel) 
or even for one set of observation (e.g. assuming that some 
parameters of BRDF are spectrally independent in any single 
pixel). 

The approach is the following. First, the Normal system 
is built under the assumption that every single parameter that 
drives the forward model is retrieved. Then, as shown in 
Appendix C, fixing several parameters equal in the inversion 
can be achieved by rather simple modification of the Normal 
system without any other changes in the inversion algorithm. 
Assume one has defined the Fisher matrix A^ and gradient 
VfldN) for deriving N parameters ai . Then if we assume that 
n + 1 parameters are equal, i.e. at =ai k (k= 1 , ..., n + 1), we 
should decrease number of the retrieved parameters and re- 
calculate Fisher matrix and gradient for new decreased set of 


the retrieved parameters. Direct realization such reduction of 
the retrieved set of parameters requires significant changes in 
the algorithm. However, analyzing structure of the equations 
we noted that the Fisher matrix and gradient for the reduced 
set of parameters can be obtained trivially from the Fisher 
matrix and gradient calculated for original extended set of 
retrieved parameters. Specifically, we can implement the re- 
trieval simply by decreasing the dimensions of the Fisher ma- 
trix and gradient by summing up n + 1 of lines and columns 
of the equal parameters. As a result the solution can be ob- 
tained by solving modified Normal system with Fisher ma- 
trix and gradient of smaller dimensions: 

(N —n) x (N —n) and (N — n) correspondingly. The de- 
tails are shown in Appendix C. It is important to note that 
the described reorganization of the Normal system is rather 
straightforward and is not time consuming. 


5 Algorithm functionality and sensitivity tests 

The algorithm is designed to retrieve a rather extended set of 
parameters describing the atmospheric aerosol and underly- 
ing surface reflectance. The list of the parameters is given 
in Table 1 . The algorithm derives basically the same group 
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of aerosol parameters in both scenarios of observation over 
ocean or land. At the same time, there is some flexibility in 
setting the representation of the aerosol size distribution. The 
size distribution can be modeled using different number of 
size bins A; : (i) large number of the size bins (A/ >10) mod- 
eled as tri-angle functions; (ii) small number of the size bins 
(A; < 10) modeled as log-normal size functions with differ- 
ent assumed width. The smaller number of bins can be set 
for reducing retrieval implementation time. The retrieval of 
large numbers of size bins can be a preferable option in sit- 
uations where observations are more sensitive to details of 
aerosol properties. For example, contribution of the atmo- 
spheric aerosol dominates the POLDER observations over 
the dark ocean surface. The modeling of ocean surface re- 
flection fully follows the approach the currently operational 
POLDER algorithm (Deuze et al., 2001; Herman el al., 2005; 
Tanre et al., 2011). The only difference assumed in the de- 
veloped algorithm is a possibility to retrieve the values of the 
Lambertian reflectance linked to the seawater reflectance at 
short wavelengths (0.44 and 0.49 pm) and to the whitecaps 
reflectance properties. The details of algorithm performance 
in this situation have not been tested yet. We plan to pro- 
vide analysis and illustrations of the retrieval algorithm per- 
formance over ocean in future publications. 

The design of the inversion algorithm allows two comple- 
mentary inversion scenarios. The conventional Single-Pixel 
Inversion is designed to retrieve all above aerosol and surface 
parameters for each single pixel of observation. The Multi- 
Pixel Retrieval implements the simultaneous inversion of 
satellite observations over a group of the neighboring pixels 
obtained during the limited time period (e.g. several weeks). 
This strategy allows applying rather diverse constraints on 
variability of every retrieved parameter in space and time. It 
also allows the assumption of time or space pixel indepen- 
dence of any single parameter or a group of the retrieved pa- 
rameters. Based on the known experience in aerosol and land 
surface retrieval developments, and on general understanding 
of limitations of the information content (checked in sensitiv- 
ity tests) the following combination of assumptions has been 
chosen for applying to POLDER/PARASOL observations. 

For aerosol (within each single pixel): 

- The aerosol size distribution is assumed smooth. 

- The spectral dependence of n(X) and k (A./) is assumed 
smooth. 

For aerosol ( between pixels): 

- The moderate spatial variability is assumed for all 
aerosol parameters with the exception of aerosol load- 
ing (driven by aerosol total concentration). 

- Significant spatial variability of aerosol loading (total 
concentration) is allowed. 


For surface reflectance (within each single pixel): 

- p 0 (X) is assumed moderately spectrally smooth. 

- k (X) and 9 (X) are assumed strongly spectrally smooth 
(nearly constant). 

- (if retrieved ) h Q (X) is assumed moderately spectrally 
smooth. 

- B (X) is assumed strongly to moderately spectrally 
smooth. 

For surface reflectance (between pixels): 

- All parameters of surface reflectance are assumed con- 
stant during ~7 days. 

- No constraints are applied at spatial variability of the 
parameters. 

Table 5 summarizes the above assumptions and Tables 3 
and 4 show the Lagrange parameters used for applying 
smoothness constraints on the size or spectral dependencies 
of the retrieved aerosol and surface properties in each single 
pixel and on the spatial and temporal variability of the any 
single retrieved parameter between different pixels. 

A series of sensitivity tests has been performed to verify 
the performance and potential of the developed algorithm. 
The tests have been designed to provide rather compact and 
conclusive illustration of capabilities and limitations of the 
algorithm to derive full set of aerosol and surface reflectance 
parameters from POLDER type satellite observations. The 
discussion below will be focused on the satellite observa- 
tions over land surfaces. The tests verifying the performance 
of the algorithm over dark surfaces are not discussed in this 
paper. Nonetheless such tests were also performed and they 
showed rather robust retrieval of all aerosol parameters in 
most tested situations. Some results of such tests have been 
shown in the paper by Kokhanovsky et al. (2010) that dis- 
cusses the outcome of the series of aerosol retrieval “blind 
test”. In the framework of this study the observations of 
aerosol over dark surfaces by different satellites were sim- 
ulated for a single chosen undisclosed aerosol model. These 
observations were distributed to different research groups in- 
volved in aerosol retrieval developments. The groups willing 
to participate in such tests have derived the aerosol proper- 
ties from the obtained synthetic observations and returned 
the obtained aerosol retrieval results to the distributor of the 
data. Once all retrieval results were collected, the assumed 
aerosol model was disclosed and compared with the collected 
retrieval results. Based on the analysis of the “blind test” 
outcome the present algorithm (mentioned in the paper by 
Kokhanovsky et al., 2010, as “LOA-2”) was rated among 
the algorithms providing the most comprehensive and ac- 
curate results. The retrieval results provided for studies by 
Kokhanovsky et al. (2010) were obtained using conventional 
pixel-by -pixel retrieval approach. 
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Table 5. Summary of the a priori assumptions applied for constraining of each element of the vector of unknowns. 


AEROSOL: 

Parameter: 

Single-Pixel 

constraints 

Inter-Pixel 

constraints 

C v 

dV( ri )/dlnr 
n(ki) 
k(ki ) 

no 

weak size smoothness 
strong spectral smoothness 
mild spectral smoothness 

t - no; X - mild; Y - mild 

t - constant during 7 days; X-Y ( horizontal ) - mild 
t - constant during 7 days; X-Y ( horizontal ) - mild 
t - constant during 7 days; X-Y ( horizontal ) - mild 

SURFACE REFLECTANCE: 

Parameter: 

Single-Pixel 

constraints 

Inter-Pixel 

constraints 

Poiki) 

K(ki) 

0(ki) 

ho(ki) 

B(ki) 

weak spectral smoothness 
strong spectral smoothness 
strong spectral smoothness 
weak spectral smoothness 
strong spectral smoothness 

t - constant during 7 days; X-Y ( horizontal ) - mild 
t - constant during 7 days; X-Y ( horizontal ) - mild 
t - constant during 7 days; X-Y ( horizontal ) - mild 
t - constant during 7 days; X-Y ( horizontal ) - mild 
t - constant during 7 days; X-Y ( horizontal ) - mild 


The sensitivity tests for retrieval over land surface were 
aimed to verify the performance of the retrieval approach un- 
der conditions maximally reproducing the real environment. 
With that purpose, the POLDER observation geometry and 
the aerosol and surface characteristics have been assumed to 
mimic closely the observations over two AERONET (Hol- 
ben et al., 1998) observation sites, Banizoumbou (Niger) and 
Mongu (Zambia). The properties of both aerosol and sur- 
face reflectance were extensively studied and discussed in 
the scientific literature. For example, the observations over 
Banizoumbou were discussed in a number of papers (Rajot 
et al., 2008, Formenti et al., 2008; Johnson et al., 2009; 
Sow et al., 2009, etc.), and the observations over Mongu 
were discussed in many studies (Dubovik et al., 2002a; Eck 
et al., 2003; Schmid et al., 2003; Haywood et al., 2003; 
Pilewskie et al., 2003; Lyapustin et al., 2006; Sinyuk et al., 
2007, etc.). Two rather distinct aerosol types dominate the 
aerosol over AERONET sites: desert dust and biomass burn- 
ing. Mongu is highly affected by biomass burning aerosol. 
The aerosol over Banizoumbou site is strongly impacted by 
desert dust outbreaks with notable seasonal contributions 
of smoke. Biomass burning is known to be dominated by 
fine absorbing spherical particles. All optical properties of 
biomass burning aerosol (extinction, single scattering albedo, 
etc.) generally have distinctly decreasing spectral depen- 
dence (e.g. Eck et al., 1999; Dubovik et al., 2002a; Reid 
et al., 2005, etc.). In contrast, desert dust is dominated by 
coarse non-spherical and weakly absorbing particles. There- 
fore, the extinction of desert dust is generally spectrally flat 
(e.g. Eck et al., 1999; Holben et al., 2001, etc.). Dust parti- 
cles generally are nearly non-absorbing in visible with some 
moderate absorption appearing at ~0.5 pm and increasing to- 
wards UV spectral range (e.g. Kaufman et al., 2001; Dubovik 
et al., 2002a; Lafon et al., 2006, etc.). Biomass fine aerosol 


particles strongly polarize the scatted light, while the polar- 
ization by coarse non-spherical desert dust particles is weak 
(e.g. Volten et al., 2001; Dubovik et al., 2006, etc.). The 
two series of tests were conducted: (i) when desert dust 
was dominating and (ii) when biomass burning was dom- 
inating. The sensitivity tests were carried out for a wide 
range of aerosol loadings. The 16 test scenarios were de- 
signed to cover the range from r (0.44 pm) = 0.01 to r (0.44 
pm) = 4. The size distributions for both biomass and desert 
dust aerosol models were adapted from original observa- 
tions over Banizoumbou and Mongu AERONET sites. The 
22 AERONET size distribution bins were used for generating 
synthetic observations. Since the size distributions were 
adapted from actual AERONET observations, they did not 
have perfect multi-modal shape. This fact makes the tests 
more challenging but more realistic. The values of complex 
refractive index for A = 0.44, 0.67, 0.87 and 1.02 pm were 
adapted from actual AERONET observations. The values 
for intermediate spectral channels k = 0.49 and 0.55 pm were 
obtained by the interpolation. At the same time it should be 
noted that values chosen for both the size distribution and 
complex refractive index are generally close to the values re- 
ported from AERONET retrieval climatology for dust and 
biomass burning by Dubovik et al. (2002a). 

The random noise at the level of 1 % for intensity and 
0.5 % for degree of linear polarization of PARASOL signal 
have been added to the simulated “synthetic” PARASOL ob- 
servations. These synthetic observations were inverted by the 
present algorithm using two retrieval scenarios: conventional 
pixel-by -pixel and multi-pixel retrieval. In the first scenario 
all synthetic observations were inverted fully independently. 
For multi-pixel retrieval, the “synthetic” PARASOL observa- 
tions were assumed to be observed during ~1 week in four 
consequent observations (with time difference of 1 to 2 days). 
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Table 6. Definition of “ the initial guess” for the vector of unknowns. 


Aerosol Properties 

Surface Reflectance 


C v = Cq (corresponding to the value of r ae r (0.44)^0.05); 

Po(kj) = 0.05 (i = l,. 

N k ) 

dV(ri)/d\nr = QA; (i = 1, ..., N r ) 

K(ki) = 0.75 6 = 1, ... 

,N k ) 

Qph — 0.7 

6>(A.;) = -0.1 6 = 1, .. 

; N k ) 

*(A. f ) = 1.4 = 

ho(^i) = Poi^i) 6 = 1 

U.,A6l) 

k(Xi) = 0.005 6 = 1,..., N k ) 

B(ki) = 0.03 6 = 1,.. 

;N k ) 





Assumed x(0.44 pm) 


Assumed x(0.44 pm) 


Assumed x(0.44 pm) 


Fig. 8. Retrieval of biomass burning optical thickness over Mongu for 16 different aerosol loadings from r (0.44 pm) = 0.01 to 4.0: (left 
panel) single-pixel retrieval with no noise added; (center panel) single-pixel retrieval with the random noise added; (right panel) multi-pixel 
retrieval with the random noise added. 


It was assumed that the four neighboring pixels were ob- 
served at each observation time. For example, Fig. 6 shows 
situation, when 9 (3 x 3, where N x = 3 and N y = 3) pixels ob- 
served during 3 consequent observations. One can imagine 
similar case when 4 (2 x 2, where N x =2 and N y = 2) pixels 
observed during 4 consequent observations. The constraints 
were applied on the inter-pixel variability of the retrieved 
aerosol and surface properties as described above in this Sec- 
tion and summarized in Tables 2-5. Specifically, the proper- 
ties of the surface reflectance were assumed constant during 
the week for the same pixels. No constraints were applied on 
either temporal or horizontal variability of aerosol concentra- 
tion and for all other aerosol parameters only weak horizontal 
variability was allowed during the same day (i.e. observed by 
the same PARASOL image). All retrievals in the tests were 
conducted using most detailed size distribution representa- 
tion: 16 logarithmically equidistant size bins covering the 
range of aerosol particle radii from 0.1 to 7 pm. The same 
initial guess was used for the unknown parameters in every 
pixel. It is described in Table 6. 

Figures 8-1 1 show the results of the sensitivity test for re- 
trieving biomass burning aerosol over Mongu. The illustra- 
tions demonstrate the retrieved normalized size distribution, 
r (0.44 pm), co Q (k) and spectral albedo of the surface calcu- 
lated using the retrieved parameters of aerosol microphysics 


and surface reflectance respectively. Each figure shows the 
retrieval results for three situations. Left part of each figure 
shows the results of the retrieval for the case with no noise 
added. It should be noted however, that some discrepancies 
were present even in this case since the size distribution is 
modeled using 22 logarithmically equidistant size bins cov- 
ering range of particle radii from 0.05 to 15 pm, while re- 
trieval uses only 16 size bins covering radii from 0.1 to 7 
pm. The central and the right parts of each figure show the 
results of the retrieval with no noise added at the level of 1 % 
for total radiances and 0.5 % for degree of linear polariza- 
tion. Both the left and central parts show the results obtained 
using conventional pixel-to-pixel retrieval approach, while 
the right part of each figure shows results obtained using the 
multi-pixel retrieval strategy. Figures 12-15 are analogous to 
Figs. 8-1 1 with the difference that they illustrate the retrieval 
results of desert dust over Banizoumbou (Niger). 

The analysis of the results shows that in situations where 
no noise is added, all aerosol characteristics are retrieved 
rather accurately. Some notable deviations are seen for sin- 
gle scattering albedo and size distribution in situations with 
very low aerosol loading when r (0.44 pm) < 0.2. The de- 
viations are particularly notable for the size distribution re- 
trievals of the coarse size particles and for single scatter- 
ing albedo of biomass burning at long wavelengths. These 
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Fig. 9. Retrieval of biomass burning size distribution over Mongu for 16 different aerosol loadings from r (0.44 pm) = 0.01 to 4.0: (left 
panel) single- pixel retrieval with no noise added; (center panel) single-pixel retrieval, with the random noise added; (right panel) multi-pixel 
retrieval, with the random noise added. 
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Fig. 10. Retrieval of biomass burning single scattering albedo over Mongu for 16 different aerosol loadings from r (0.44 pm) = 0.01 to 
4.0: (left panel) single-pixel retrieval with no noise added; (center panel) single-pixel retrieval with the random noise added; (right panel) 
multi-pixel retrieval with the random noise added. 
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Fig. 11. Retrieval of surface albedo over Mongu for 16 different aerosol loadings from r (0.44 pm) = 0.01 to 4.0: (left panel) single-pixel 
retrieval with no noise added; (center panel) single-pixel retrieval with the random noise added; (right panel) multi-pixel retrieval with the 
random noise added. 
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Fig. 12. Retrieval of desert dust optical thickness over Banizoumbou for 16 different aerosol loadings from r (0.44 pm) = 0.01 to 4.0: (left 
panel) single-pixel retrieval with no noise added; (center panel) single-pixel retrieval with the random noise added; (right panel) multi-pixel 
retrieval with the random noise added. 



Fig. 13. Retrieval of desert dust size distribution over Banizoumbou for 16 different aerosol loadings from r (0.44 pm) = 0.01 to 4.0: (left 
panel) single pixel retrieval with no noise added; (center panel) single-pixel retrieval, with the random noise added; (right panel) multi-pixel 
retrieval with the random noise added. 


deviations can be explained by the fact that at such low 
aerosol loadings the aerosol contribution into radiation ob- 
served by satellite is negligible compare to the contribution 
of land surface reflectance. Therefore even minor perturba- 
tions of the observation may significantly affect the retrieval 
results. In contrast, in situations with very high aerosol load- 
ing r (0.44 pm) > ~ 2-2.5, the contribution of aerosol dom- 
inates in the observed reflected radiation and the reflectance 
properties of the underlying land surface become nearly in- 
visible for satellite. Correspondingly, the retrievals of surface 
reflectance become unstable, as it is seen for retrieved sur- 
face reflectance values at short wavelengths. All these ten- 
dencies become pronounced once the random noise is added. 
As one can see from central parts of Figs. 8-15, the size 
distribution and single scattering albedo retrieved by con- 
ventional pixel-by-pixel approach deteriorate in most of the 
cases including the situations with moderate and even high 
aerosol loading with r (0.44 pm) reaching 0.4 and even 0.8 
(for single scattering albedo in particular). In addition one 


can see the significant complications in the retrieval of shape 
of aerosol size distribution, in particular for larger particle 
sizes. The size distribution retrieved for desert dust (see 
Fig. 13) becomes much wider and the maximum shifts to- 
wards smaller sizes. In retrievals of biomass burning size dis- 
tributions (Fig. 9) there are also significant complications for 
the large radii. This difficulty can be explained by the well- 
known fact (e.g. see Bohren and Huffman, 1983) of strongly 
decreasing scattering efficiency (per unit of particle vol- 
ume, i.e. k ext (X,r)/v(r)) for particles with sizes much larger 
than the wavelength of observations. Therefore retrieving 
accurate shape of size distribution for radii > ~3 pm from 
POLDER/PARASOL observations appears to be very diffi- 
cult. The retrieval of aerosol optical thickness r (0.44 pm) 
seems to be rather reliable in situations when random noise is 
added with the exception of cases of very high aerosol load- 
ing r (0.44 pm) > ~3 where the retrieval errors reach the val- 
ues of 0.5 or even larger (see Fig. 8). These high errors can 
be explained by the fact that in such situations the reflected 
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Fig. 14. Retrieval of desert dust single scattering albedo over Banizoumbou for 16 different aerosol loadings from r (0.44 pm) = 0.01 to 
4.0: (left panel) single-pixel retrieval with no noise added; (center panel) single-pixel retrieval with the random noise added; (right panel) 
multi-pixel retrieval with the random noise added. 
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Fig. 15. Retrieval of surface albedo over Banizoumbou for 16 different aerosol loadings from r (0.44 pm) = 0.01 to 4.0: (left panel) single- 
pixel retrieval with no noise added; (center panel) single-pixel retrieval with the random noise added; (right panel) multi-pixel retrieval with 
the random noise added. 


radiances are dominated by the multiple scattering of very 
high orders and some of fine agular features of reflected light 
distribution are smoothed out. Correspondingly, deriving 
detailed aerosol and surface properties becomes more dif- 
ficult. The analysis of the right parts of Figs. 8-15 shows 
that using multi-pixel approach significantly improves the 
retrievals of all aerosol and surface parameters. This ten- 
dency is not surprising because added constraints allow prop- 
agation and consolidation of useful information from differ- 
ent observational situations. For example, in the situations 
with low aerosol loading the satellite observes mainly sur- 
face reflectance properties. Correspondingly, once the con- 
straints limiting time variability of the surface reflectance are 
applied, this information is supplied into the interpretations 
of observations corresponding to moderate and high aerosol 
loading over the same pixel. Similarly, the constraints of 
horizontal variability of aerosol properties help to improve 
the retrieval of aerosol by benefiting from observations of 


the same or similar aerosol properties over several pixels 
with somewhat different conditions of observations (geom- 
etry, surface reflectance, aerosol loading). 

It should be noted that such retrieved parameters as aerosol 
mean height, fraction of spherical particles, detailed param- 
eters of BRDF and BPDF are not shown in Figs. 8-15. All 
these parameters have been included in the tests. The re- 
sults are shown only for most significant aerosol and surface 
parameters commonly discussed in remote sensing retrieval 
analysis. In a future study, it is planned to implement the 
comprehensive sensitivity analysis of the retrieval approach 
suggested here. The sensitivity studies shown in the present 
paper are aimed to provide insightful but preliminary outlook 
at the expected performance of the retrieval. Nonetheless, it 
is possible to state here that the tests have shown reasonably 
robust retrieval of all sought parameters. Figures 11 and 15 
show the retrieval of surface albedo. The retrievals of the 
detailed BRDF and BPDF parameters generally demonstrate 
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Fig. 16. Statistics of the retrieval errors of biomass burning optical thickness over Mongu in the series of the numerical tests for 100 different 
realizations of random noise added (multi-pixel retrieval): (left panel) 0.44 pm; (right panel) 1.02 pm. 


the same trends. In pixel-by -pixel retrieval scenario all pa- 
rameters were retrieved rather accurately in situations with 
low and moderate aerosol loading. Using multi-pixel re- 
trieval helped to achieve stable retrieval of all surface re- 
flectance parameters in all situations including the cases 
with high aerosol loading. In contrast, the retrieval of the 
mean altitude of the aerosol layer h a was retrieved robustly 
when aerosol loading was moderate or high. When pixel- 
by-pixel retrieval scenario was used, h a was retrieved with 
accuracy better than 1 km for r (0.44 pm) > ~0.5 and bet- 
ter than 0.5 km for r (0.44 pm) > ~1.0. Once the multi-pixel 
scenario was applied, h a was retrieved with accuracy better 
than 0.5 km for r (0.44 pm) > ~0.5 and better than 0.3 km 
for r (0.44 pm) > ~ 1 .0. The retrieval of the spherical particle 
fraction was rather successful for desert dust aerosol, where 
the coarse mode is dominant. When pixel-by -pixel retrieval 
scenario was used, C sp h was retrieved with the accuracy of ~ 
50 % in situation with r (0.44pm) >~1.0. In multi-pixel re- 
trieval scenario it was retrieved with the accuracy of ~30 % 
in situations with r (0.44 pm) > ~0.2 and better than 10% 
for r (0.44 pm) > ~1.0. 

Finally, Figs. 16-19 summarize the performance of the re- 
trieval in an extended series of numerical tests with added 
random noise. The series were composed by the tests anal- 
ogous to those illustrated in Figs. 8-15, with the difference 
that each test was repeated 100 times with different realiza- 
tions of the generated noise. The random noise was added at 
the level of standard deviation cr = 1 % for total radiances and 
cr =0.5 % for degree of linear polarization. Thus Figs, lb- 
19 show the results of 3200 inversions = 100 x 16x2 (100 

- number of tests with different realizatios of modeled 
rundom noise; 16 - aerosol loading covering r (0.44 pm) 
from 0.01 to 4.0; 2- number of observational sites: Ban- 
izoumbou and Mongu) and the differences: (r (retrieved) 

- r (“true”))! r (“true”) and co 0 (retrieved) - <u w (“true”) for 


X = 0.44 and 1.02 pm. All retrievals were conducted using 
most robust multi-pixel retrieval scenario. Table 7 provides 
a brief summary of the test outcome. It provides the stan- 
dard deviations cr a t /t (A.) of the relative errors calculated us- 
ing the retrievals obtained with 100 different realizations of 
added random errors as follows: 


^Ar/i 


and 


\ 


- T 

N 4 -^= 1 , 


^retrieved 


^true 


rtrue 


x 100 (%) (35a) 


; _ | m f^O. retrieved ^0,true^ • (35b) 


These estimates can be considered as reference values for the 
expected retrieval errors. It should be noted that the real er- 
rors in actual POLDER/PARASOL observations (Fougnie et 
al., 2007) are likely to have 2 to 3 times higher levels than the 
errors modeled. Therefore, the actual retrieval errors should 
be expected to be higher than the numbers given in the Ta- 
ble 7 by a factor of 2 or 3. In any case, these estimates are to 
be verified in a more comprehensive sensitivity analysis. 


6 Algorithm applications to real POLDER/PARASOL 
observations 

As a final stage of this study the developed algorithm 
was applied to actual observations by POLDER/PARASOL. 
For consistency with the sensitivity test, the algorithm was 
used to process the full 2009 year of POLDER/PARASOL 
data over Banizoumbou (Niger) and Mongu (Zambia) 
AERONET sites. Specifically, the algorithm was applied 
to 4 POLDER/PARASOL pixels surrounding the exact lo- 
cations of the sites. The multi-pixel inversion scenario was 
employed in the retrieval. Since the current studies did not 
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Fig. 17. Statistics of the retrieval errors of biomass burning single scattering albedo over Mongu in the series of the numerical tests for 
100 different realizations of random noise added (multi-pixel retrieval): (left panel) 0.44 pm; (right panel) 1.02 pm. 
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Fig. 18. Statistics of the retrieval errors of desert dust optical thickness over Banizoumbou in the series of the numerical tests for 100 different 
realizations of random noise added (multi-pixel retrieval): (left panel) 0.44 pm; (right panel) 1.02 pm. 


address the cloud-screening of the satellite data, the algo- 
rithm was applied to PARASOL data identified as cloud-free 
by the current operational POLDER algorithm (Deuze et al., 
2001; Herman el al., 2005; Tanre et al., 2011). In addition, 
the retrieval outliers with the fitting residual higher than 5 % 
were eliminated from the final results. 

It should be noted that the algorithm has been developed 
here for utilization of complete set of POLDER/PARASOL 
observations as shown in Table 1. However, the analysis 
of PARASOL in-flight calibration and performance studies 
performed by Fougnie et al. (2007) indicate that PARASOL 
0.443 pm band does not meet to the mission requirements 
due to unidentified stray light problem and, therefore, it is 
not recommended for use. Nonetheless for consistency with 
the sensitivity tests performed in Section 5, we have de- 
cided to use PARASOL measurements at 0.443 pm in these 


illustrative applications of the algorithm. Such a decision 
was also made with the idea to benefit from unique informa- 
tion content of the observations at the shortest wavelength. 
Indeed, the optical thickness of aerosol is generally higher at 
0.443 pm compare to other PARASOL spectral bands, while 
the land surface reflectance is generally lower at 0.443 pm 
compare to other bands. In order to minimize the possible 
effect of the higher uncertainty at 0.443 pm spectral band, we 
have introduced the coefficient for correcting observations at 
this wavelength. Specifically, the data were inverted using 
a set of different correction coefficients that could increase 
or decrease the measured radiances at 0.443 pm up to 10 %. 
The coefficient corresponding to the best fit of corrected ob- 
servations by the theoretical model was used in the illustra- 
tions below. Such strategy for correction of PARASOL ob- 
servation is in agreement with evaluation of possible effect of 
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Table 7. The summary of the sensitivity tests for aerosol retrievals from simulations mimicking POLDER observations over Mongu (Zambia) 
and over Banizoumbou (Niger) with the added random noise at the level of standard deviation a = 1% for total radiances and a = 0.5% for 
the degree of linear polarization. 
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Fig. 19. Statistics of the retrieval errors of desert dust single scattering albedo over Banizoumbou in the series of the numerical tests for 
100 different realizations of random noise added (multi-pixel retrieval): (left panel) 0.44 pm; (right panel) 1.02 pm. 


stray light contamination (B. Fougnie, personal communica- 
tion, 2011). In future studies we plan to investigate compre- 
hensively the propagation of the PARASOL 0.443 pm band 
errors to the aerosol retrieval results. We hope to identify 
the best strategy for addressing this uncertainty in the oper- 
ational application of the current algorithm of using the cor- 
rected 0.443 pm data (e.g. as it is done here) or eliminating 
them completely from the aerosol retrieval. 

The r (0.44 pm) and co Q (0.44 pm) calculated using the set 
of retrieved aerosol parameters (Table 2) were compared with 
available AERONET data. The comparisons showed rather 
robust performance of the algorithm. The retrieved values 
of aerosol optical thickness closely followed the AERONET 
observations with a correlation coefficient of ~0.9 for Ban- 
izoumbou and ~0.87 for Mongu. The values of aerosol 
single scattering albedo also agreed reasonably well with 
AERONET data for observations with high aerosol loading. 
The differences between the values of co 0 (0.44 pm) obtained 


from PARASOL and from AERONET did not exceed 0.03- 
0.05 for the cases when r (0.44 pm) > ~0.5. It also should 
be noted that the retrieval with different correction coeffi- 
cients did not exhibit dramatic changes in the retrieval. For 
example, the differences in retrieval of r (0.44 pm) and co 0 
(0.44 pm) due to changes in the correction coefficient had 
generally smaller magnitudes than mean differences between 
PARASOL and AERONET results. 

Figures 20-21 illustrate the comparison of POLDER/ 
PARASOL retrievals with the AERONET observations dur- 
ing 2 months over Banizoumbou and Mongu. These particu- 
lar periods were chosen for the illustrations because they cor- 
responded to the most complete and continuous AERONET 
data records during the seasons with the presence of high 
aerosol loadings. As can be seen from the illustrations, the 
aerosol loading trends retrieved from PARASOL agree well 
with the AERONET observations. The values of aerosol sin- 
gle scattering albedo are also in reasonable agreement with 
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Fig. 20. The comparison of r (0.44 pm) and coq (0.44 pm) retrieved from POLDER/PARASOL during August-September 2009 over Mongu 
with the corresponding values provided by AERONET. 
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Fig. 21. The comparison of r (0.44 pm) and coq (0.44 pm) retrieved from POLDER/PARASOL during January-February 2009 over Bani- 
zoumbou with the corresponding values provided by AERONET. 


AERONET values when r (0.44 pm) ~0.5 or larger. At the 
same time, one can identify some notable differences be- 
tween AERONET data and PARASOL retrievals in some 
single points. Nonetheless, our analysis showed that usually 
these points correspond to days when AERONET data indi- 
cate the partial cloudiness during the day as identified by the 
Smirnov et al. (2000) cloud-screening procedure. Therefore, 
such outliers can probably be explained by some sky inhomo- 
geneities in one or all PARASOL observed pixels that cover 
an area of 12 x 12 km around the AERONET site. It should 
be also noted that comparisons of AERONET observations 
and PARASOL retrieval results showed that the retrievals 
tend to underestimate the spectral dependence of retrieved 
aerosol properties. The desert dust retrievals tend to show 


more pronounced spectral dependence of aerosol properties, 
while the aerosol properties retrieved for biomass burning 
tend to have smaller spectral dependence than observed from 
AERONET. As follows from the sensitivity tests and our gen- 
eral understanding, this is likely due to very low sensitivity 
of satellite observations to the shape of the aerosol size dis- 
tribution for large particles. Indeed, the contribution of very 
large particles (r >3 pm) into radiation reflected to space is 
significantly smaller than from fine particles, while the spec- 
tral variability of aerosol properties is very sensitive to the 
presence of such large particles. Therefore, deriving fully 
adequate spectral dependence of aerosol optical properties 
from satellite observations appears to be a harder task than 
deriving aerosol loading especially over bright land surfaces 
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considered in present test. Apart from this relatively minor 
issue, the present analysis showed very encouraging perfor- 
mance of the developed algorithm over bright land surfaces, 
i.e. in conditions that are considered traditionally the most 
difficult for retrieval of aerosol from satellites. This result is 
particularly encouraging because the algorithm is designed to 
provide rather extensive set of the retrieved parameters pro- 
viding detailed characterization of the properties of aerosol 
and the underlying surface. Both the results of numerical 
sensitivity tests and the obtained results of actual PARA- 
SOL data inversion suggest that the processing the PARA- 
SOL with developed algorithm can provide global products 
of total aerosol optical thickness and single scattering albedo 
rather consistent with the accuracy requirements formulated 
by Mishchenko et al. (2004, 2007). However, some limita- 
tions in accuracy and scope of the aerosol information de- 
rived by presented algorithm should be noted. For example, 
the present approach does not discriminate between optical 
properties of aerosol fine and coarse modes (see Sect. 3.3). 
Such retrieval assumption is consistent with the resutls of 
limited tests conducted in scope of this study for POLDER 
type imager. In addition, analyzing the standard deviations 
cfAr/r (k) resulted of retrieval tests with synthetic PARASOL 
observations perturbed of random noise and shown in Table 7 
one can see that errors in the retrieved total aerosol optical 
thickness r (k) (especially for high values of r) could ex- 
ceed the value of 0.04 suggested by Mishchenko et al. (2004, 
2007) as APS retrieval uncertainty expected over land for 
optical thickness of both fine and coarse modes of aerosol. 
At the same time, the results presented here are preliminary. 
Further testing, verification and tuning of the presented algo- 
rithm are planned. In addition, the efforts are also planned 
for the acceleration of the developed algorithm. In the cur- 
rent state of the algorithm the implementation of the retrieval 
for a single PARASOL pixel (5.3 x 6.2 km) required in av- 
erage of ~10s of computer time. This time should be de- 
creased in 50 to 100 times, in order to reprocess existing data 
archive of PARASOL observations in reasonable time frame- 
work. Such acceleration of the algorithm is expected to be 
achieved by reengineering and optimizing the algorithm, for 
example, using parallel programming and implementing the 
retrievals at several computers in parallel. 


7 Conclusions 

The paper has discussed in detail a concept for a new state- 
of-the-art algorithm developed for deriving detailed proper- 
ties of atmospheric aerosol from satellite observations. The 
proposed retrieval does not use precalculated look-up ta- 
bles commonly utilized in satellite retrievals for fitting ob- 
servations. Instead, a more general approach is applied that 
searches in continuous space for the solutions and optimizes 
the statistical properties of the obtained retrieval. Such op- 
timization can be achieved by adjusting the structure of the 


Table 8. Abbreviations. 


AERONET 

APS 

AVHRR 

BRDF 

BPDF 

LSM 

MERIS 

MISR 

MODIS 

MSG 

PARASOL 

POLDER 

RSP 

SEVIRI 

TOMS 


AErosol RObothic NETwork 
Aerosol Polarimetry Sensor 
Advanced Very High Resolution Radiometer 
Bidirectional Reflectance Distribution Function 
Bidirectional Polarization Distribution Function 
Least Square Method 

Medium Resolution Imaging Spectrometer Instrument 
Multiangle Imaging SpectroRadiometer 
Moderate Resolution Imaging Spectroradiometer 
Meteosat Second Generation 

Polarization and Anisotropy of Reflectances for Atmospheric 
Science coupled with Observations from a Lidar 
POLarization and Directionality of the Earth’s Reflectances 
Research Scanning airborne Polarimeter 
Spinning Enhanced Visible and InfraRed Imager 
Total Ozone Mapping Spectrometer 


deviations in an effort to fit observations by theoretical model 
under conditions where the amount of observations exceeds 
the number of retrieved parameters. The set of observations 
provided by modern enhanced spectral multi-viewing spec- 
tral polarimeters allows applying such optimization. The 
algorithm described in this paper was adapted for applying 
to observations of POLDER/PARASOL imager that regis- 
ters reflected atmospheric radiation at six wavelengths in up 
to 16 directions. The algorithm fits total radiance and lin- 
ear pollarization observed in all directions in all available 
spectral channels using generalized multi-term Least Square 
type numerical inversion formulation. That formulation al- 
lows fitting several sets of both observations and a priori data. 
The concept complementarily unites advantages of a variety 
of practical inversion approaches, such as Phillips -Tikhonov - 
Twomey constrained inversion, Kalman filter, Newton-Gauss 
and Levenberg-Marquardt iterations, etc. This methodology 
has resulted from a multi-year effort at developing inver- 
sion algorithms for retrieving comprehensive aerosol prop- 
erties from AERONET ground-based observations. The al- 
gorithm is driven by a large number of unknowns and de- 
signed as a retrieval of extended set of parameters affecting 
measured radiation. For example, over land the algorithm is 
set to retrieve parameters of the land surface reflectance to- 
gether with detailed information about aerosol sizes, shape, 
absorption and composition (refractive index) and aerosol 
layer elevation. 

In addition, the algorithm is developed as simultaneous in- 
version of a large group of pixels within one or several im- 
ages. Such, multi-pixel retrieval regime takes advantage of 
known limitations on spatial and temporal variability in both 
aerosol and surface properties. Specifically the pixel-to-pixel 
and/or day-to-day variations of the retrieved parameters are 
enforced to be smooth by additional appropriately set a priori 
constraints. This concept is aimed to achieve higher consis- 
tency in satellite retrievals because in such an approach the 
solution over each single pixel is benefiting from information 
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contained in coincident observations over neighboring pix- 
els, as well as from information about surface reflectance 
(over land) obtained in preceding and consequent observa- 
tions over the same pixels. The paper provided detailed de- 
scription of full set of formulations necessary for realizing 
this concept. 

The performance of the developed algorithm has been 
demonstrated by application to both synthetically generated 
and real POLDER/PARASOL observations. First, a series 
of sensitivity tests was conducted by applying the algorithm 
to the synthetic POLDER/PARASOL observations over veg- 
etated and desert surfaces. The simulations were designed 
to mimic satellite observations over well-studied AERONET 
network sites in Mongu (Zambia) and Banizoumbou (Niger). 
The synthetic POLDER/PARASOL signals were perturbed 
by random noise prior to applying the retrieval algorithm. 
Both the conventional pixel-by -pixel and newly suggested 
multi-pixel retrieval approaches were tested. The results of 
the tests showed that the complete set of aerosol parameters 
can be robustly derived with acceptable accuracy in both situ- 
ations over both vegetated and desert surfaces. The summary 
of the error analysis is provided. In addition, the algorithm 
was applied to one year of PARASOL observations over both 
Mongu and Banizoumbou AERONET sites. The compari- 
son of the derived aerosol properties with available observa- 
tions by AERONET ground-based sun/sky -radiometers indi- 
cated encouraging consistency of PARASOL derived optical 
thickness and single scattering albedo with those obtained by 
AERONET. At the same time, the presented tests and anal- 
ysis of the retrieval from actual PARASOL observation had 
somewhat limited character and were aimed to provide an in- 
troduction and some limited illustration of the proposed re- 
trieval algorithm. More comprehensive studies for testing 
and tuning the developed algorithm are planned in future ef- 
forts. Such important aspects of algorithm implementation 
for operational processing as cloud- screening and retrieval 
time requirements are to be addressed in follow-on studies. 

It should be noted that the research efforts described in this 
paper considerably relied on the accumulated experience and 
many aspects of the retrieval, as well as actual computer tools 
that were inherited from precedent efforts on development 
of currently operating AERONET retrieval and PARASOL 
aerosol retrieval algorithms. 


Appendix A 

Matrices of multi-pixel smoothness constraints 

In order to determine the index of pixels we will follow first 
the changes of x,- - coordinates (i = 1, ..., N x ) than yi - coor- 
dinates (/ = 1, ..., Ny) and 0 - time coordinates (i = 1, ..., N t ). 
For the simple case when N x = N y = N t = 2, the total vector 


of unknowns a is composed of vectors of each i - th pixel 
as: 


( a (xi; yi; ti)^ 
a (x 2 ; yi; h) 
a (xi; y 2 ; *0 
a (x 2 ; y2; *i) 
a (xi; yi; t 2 ) 
a (x 2 ; yw t 2 ) 
a (xi; y 2 ; t 2 ) 
\a (x 2 ; y 2 ; t 2 ) / 


(Al) 


The spatial and temporal smoothness constraints are applied 
separately along corresponding coordinates, i.e. constraints 
of variability over x; - coordinates are applied only to the 
vectors with the same values of yi and // ; constraints of vari- 
ability over yi - coordinates are applied only to the vectors 
with the same values of x/ and ti and constraints of vari- 
ability over ti - coordinates are applied only to the vectors 
with the same values of x; and y;. Therefore, if Ax =x;+i — 
Xi = const, Ay = yt+i— yi = const and At = ti+\ — ti = const, 
for the vector a defined by Eq. (Al) the corresponding ma- 
trices of first differences are defined as: 


/1 -10 0 0 0 0 0 \ 

0 0 1 -1 0 0 0 0 

0 0 0 0 1 -1 0 0 

\0 0 0 0 0 0 1 — 1 / 


/I 0 -1 0000 0 \ 
0 1 0 -1 00 0 0 
00 0 0 1 0 -1 0 
\00 0 0 0 1 0 -1/ 


;(A2) 


/I 0 -1 0000 0 \ 

0 1 0 -1 00 0 0 
00 0 0 1 0 -1 0 
\00 0 0 0 1 0 -1/ 


Correspondingly the terms of matrix filter can be written in 
the form of the following array matrices: 



sfs, 


(( Idn 0 

V 0 Idn 
fld 21 0 

V V 0 ^ d 21 



where the smoothness matrix D is defined as: 



where Si is the matrix of first differences. The matrix I </. . is 
a diagonal matrix of dimension N d x N a (A a is dimension of 
corresponding vector a{) with the all elements on the diago- 
nal equal to corresponding element dij of matrix D, i.e.: 

Id ij = dij I(/V a x 7V a ) • (A5) 

The definition of Winter can be generalized for situations with 
a larger number of pixels. For example, for the case similar 
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to above but for N x = N y = N t = 3, Eq. (A3) are transformed 
to the following: 


sjs 
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I Id 23 0 0 \ 
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where the diagonal matrices I d tj are defined according to 
Eq. (A5) with the differences that the smoothness matrix D 
is generated as the base of the matrix of first Si or second 
differences S 2 for N x = N y = N t = 3 as: 


D = Sf Si 


/ 1 -1 0 
- 12-1 
\ 0 -1 1 


or 


(A7) 


D = Sf S 2 = 



generally significantly smaller than all the elements of Sinter 
and thus Sj nter is essentially a sparse matrix. 

It should be noted that all of the above formulations 
are given for the idealized situation where Av; = jc;+i — 
Xi - const, Ay; = y; +1 — y; = const, At; = f ;+ 1 — t; = const and 
N x = N y = N t . In reality, none of these conditions are as- 
sured and the algorithm should be able to handle the more 
general situation where A jt; / Ay; / At; / const as well as 
N x # N y ^ N t . In such general situations, the above equa- 
tions describing the matrix Sinter loose some of their simplic- 
ity and transparency. Nonetheless, the general sparse struc- 
ture of the matrix Sinter is always conserved. Also, even the 
equations loose some transparency, the realization of inter- 
pixel smoothness constraints Sj nte r is always rather simple 
on the algorithmic level and the most general case has been 
realized in the algorithm described here. 


Appendix B 

Constraining multiple-pixel retrieval by available 
information in the neighborhood of the inverted pixel 
group 

If the values of a* hr parameters a in the pixels neighboring 
the inverted pixel group are known, then the basic system 
Eq. (32) can be complemented by the additional equations: 

/* = /(«) + A / 

= ^x,y,t Cl A (A a ) 

^*,edge = ^x,edge a + A (A x a) ? (Bl) 

s *,edge = ^edge a + A (A-y a ) 

. Sledge = Sledge Cl + A (A; a) 

where S X; edge> Sledge* Sledge are matrices containing the co- 
efficients defining the m-th finite differences of parameters 
a describing the properties of the inverted pixel group with 
fl* br - These matrices can be trivially derived from S m ma- 
trices (cf. Eq. 25b). For example, in the simple case where 
N x = 3, N y = 1 and N t = 1, one can write: 


The two examples given above show a clear pattern in defin- 
ing the matrix filter for N x =N y = N t =2 and 3 that can eas- 
ily be generalized to the situation with higher dimensions. It 
can be observed from this pattern that filter always retains 
the sparse array structure. Specifically, the dimension of the 
matrix Sinter is: 

(jV a x N x x Ny x Nt) x (jv a x N x x N y x N , ), (A8) 

while the number of non-zero elements does not exceed the 
following value: 

(N a x N x x N y x N t ) x (2 x m + 1), (A9) 


( ^before 
^ after 


( a (xi; yi; 0)\ ( < 2 * \ 

a (* 2 ; y\\ h) I I a\ I 

a (xy, y\\ t\) ) \a% ) 

( a (* 4 ; yi; *i)\ ( aA \ 

a (x 5 ; yi; t\) I = <25 I 

a (x 6 ; yi; t\) J \a 6 ) 

( a (x 7 ; yt; t\)\ (a%\ 

a (* 8 ; yi; h) I I J 

a (x 9 ; yi; ti) J \a$ J 


(B2) 


If Eq. (Bl) uses the second differences, continuity with 
^before can expressed using the following formulations: 


where m is the order of the differences/derivatives used for 0 ^ = — 2 a\ + <24 + Ai 

inter-pixel smoothing. Obviously, the value of Eq. (A9) is = aj — 2 <24 + <25 + A 2 ' 


(B3) 
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In matrix form this equation can be expressed as: 

,B4 > 

Then it can be transformed as 
and one can write: 


s x ,b — ^ x,b d A (A x a), (B6) 

where 

s *x,b = - ( 0 o ~\ ) “before and = ( \ j q ) ' ( B7 > 

Analogously, the continuity with «* fter can be expressed as 
follows: 


1 0* = <25 - 2 <26 + <2y + A3 
O4 = <26 — 2 <2y + <2g + A4 

In matrix form this equation can be written as: 


(B8) 


( 


0; 

0: 


) 


/0 1 -2 1 00\ / a \ /A 3 \ 

\0 0 1 -2 10^< ft J + \ A4 / 


Then it can be transformed as 


(B9) 


/0*\ _ (0 1 -2) (l 00\ * 

\oo 1 ) “ + \ -2 1 0 ) “ after 



(BIO) 


and one can write 


< a = S x ,a« + A (A x «), (B 11) 

where 


Based on the additional equations shown in Eq. (Bl), the 
constraining of the multi-pixel system solution by additional 
conditions of continuity of aerosol and surface properties 
with values in the neighborhood of the inverted pixel group 
can be implemented by adding extra terms in the Eq. (20). 
The quadratic form (a p ), its gradient V4> (a p ) and Fisher 
Matrix A p formulated for inversion of combined multi-pixel 
Eq. (26) should be modified as follows: 



+ ^edge (# P ) > 


where 

2 ^edge (« P ) = (Sedge ~ 5 edge) ($edge flP — s edge)^B16) 


s 


edge — 


/ 1/2 „ \ 
Yx ^ a, edge 1 


/ 1/2 \ 
Yx Sx.edge \ 

1/2 „ 

Yy ‘Jy.edge 

ands edge = 

1/2 

Yy 5 y, edge 1 

1/2 c 

\ Yt a /,e dge / 


1/2 I 

\ Yt Sledge / 


.(B17) 


The Fisher Matrix A p is modified as: 


A p — 


/Ai, p 0 . 
0 A 2 , p . 

.. 0 \ 

.. 0 

T Winter T ^edge 



^ 0 0 . 

•• Atv )P j 



, (B18) 


where £2 e dge is equal to the edge smoothing inter-pixel matrix 
determined as 


S*,a = (-2 1o) and <* = -(oO l 2 ) “after' < B12 > 

Finally, using Eqs. (30) and (34), we can define s* edge and 
S x ,edge used i n Eq. (Bl) as follows: 

* -(!;;). ("if) 

Thus, the derivations shown by Eqs. (B 1)-(B 11) demonstrate 
the principle of defining the constraints that enforces the con- 
tinuity of the retrieved properties with those obtained for the 
pixels neighboring the inverted pixel group. Analogously the 
constraints can be included for continuity over coordinates y 
and t. In the general case of a large inverted pixel group one 
can use rather large numbers of equations similar to those 
of Eqs. (B3) and (B8). The maximum number of such edge 
continuity equations can be estimated as: 

N a x (N x x N y + N x x N t + N y x N t ) x (2 x m), (B14) 

where m is the order of the differences/derivatives used for 
inter-pixel smoothing. 


and the gradient (« p ) is given by 


V 4/ (a?) 


' v 4-1 (oj) 

V 4-2 (4) 


Winter — ^ 4“ V 4- e dge ( — ^ ) 


, (B20) 


where 

^ ^edge (# P ) — Sj dge ^S ec [ge Cl P — -S^ge) . (B21) 


It should be noted that the illustrative derivations shown 
in Eqs. (B3)-(B12) are given for the idealized situation 
where A Xi = %*+ 1 — xt = const and N y = N t = 1. In reality, 
the edge continuity can be applied over all three coordinates 
x,y,t and the algorithm should be enabled to handle the most 
general situation when A x; / A y ? - / A ti ^ const as well as 
N x / N y / N t . In such general situations, the above equa- 
tions describing the matrix £2 e dge loose their simplicity and 
transparency. Nonetheless, the general sparse structure of 
the matrix ft e dge is always conserved. Also, even the equa- 
tions loose some transparency, the realization of inter-pixel 
smoothness edge constraints S e dge is always rather simple on 
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the algorithmic level and it takes negligible computer time to 
implement. 

The values of the Lagrange parameters y x , y y and y t used 
for applying the “edge” inter-pixel constraints are exactly 
the same as those used in Sect. 4.2.1 describing inter-pixel 
smoothness constraints application in the simultaneous in- 
version of a group of pixels. 


Appendix C 

Assumption of inter-pixels constant parameters in 
the retrieval 

The general idea of this technique of imposing an extra as- 
sumption of equal parameters can be demonstrated using the 
following simple illustration. Assuming, Eq. (17) is repre- 
sented by the following simple system 

U « = /* + A /, (Cl) 

corresponding Normal system Aa = V4* ( a ) is defined as 

U r Ufl = U 7 '/*. (C2) 

If two parameters a\ and <22 are retrieved from one single 
observation /*, then one can write: 

u = (Ml M 2 ); f* = (/*); a = (C3) 

For such situation the Fisher matrix A and gradient V4> (a) 
are 

A® = U r U = ( “i M1 2 M2 ); (C4) 

\ U2 u\ U 2 ) 

vvpo = u r r = (Zf*)- 

If only it is assumed that a\ = <22 and only one parameter re- 
trieved, then the definitions given by Eq. (C3) should be re- 
placed by: 

U=(« 1 +B 2 );/* = (/ 1 *);« = (ai). (C5) 

Correspondingly, the Fisher matrix A and gradient V'P(fl) 
become 

A (1) = U r U = («! + m 2 ) 2 ; (C6) 

V * (1 > = U T f* = (mi + m 2 ) /*. 

Comparing Eqs. (C4) and (C5) one can note the following 
trivial relation between A^\ and A^ 2 \ 

A (1) = (A® + A®) + (A® + A®); (C7) 

V 4/ (1) = V + V 


These relationships can be easily generalized. Let us as- 
sume that we have defined the Fisher matrix A^ and gra- 
dient for deriving N different parameters < 2 /. Then 

it was decided (or discovered) that n + 1 of those initially 
different N parameters represent the same single parameter, 
i.e. ai =ai k ( k = 1, ..., n + 1). In this situation, the new Fisher 
matrix A^ N ~ n ^ and gradient V^( N ~ n ) can be obtained as fol- 
lows: 

A m,N-n = A N-nin = E A S “ non-diagonal elements, (C8) 

k=l,..„n 

^N-nh-n= E ( E A hi) “diagonal elements, (C9) 

v >1/^-") = e v ■ (cio) 

k=\,..,,n 

The resulting matrix A (N ~ n ^ has reduced dimension ( N — 
n) x (N — n). 

Thus, Eqs. (C8)-(C10) allow significant flexibility in de- 
signing the retrieval. The option is fully implemented in the 
algorithm presented here. 
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